Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 3972) +++ trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 3973) @@ -20,4 +20,24 @@ 2004/05/03: Javier Rico + * library/MSrcTranslate.[h,cc] + - added + + * library/Makefile, library/IFAELinkDef.h + - add MSrcTranslate + + * library/MSrcRotate.[h.cc] + - adapted to changes in MObservatory and MAstro + + * programs/falseSource.cc programs/falsesource.datacard + programs/srcPos.cc programs/srcpos.datacard + - added + + * programs/Makefile + - add falseSource and srcPos + + * macros/signal.C macros/signalPoint.C + - added + + 2004/05/03: Javier Rico * makeHillas.cc - put back Jobs for pedestal and calibration parts @@ -29,5 +49,5 @@ * library,programs,macros,*: - new directories added - - directory macros contains macros (.CC) and scripts (.sh) + - directory macros contains macros (.C) and scripts (.sh) - directory programs contains root compilables and associated [input] files Index: trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h (revision 3972) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h (revision 3973) @@ -9,4 +9,6 @@ #pragma link C++ class MPSFFitCalc+; #pragma link C++ class MSrcPosFromFile+; +#pragma link C++ class MSrcRotate+; +#pragma link C++ class MSrcTranslate+; #endif Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc (revision 3972) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc (revision 3973) @@ -48,7 +48,6 @@ #include "MDCA.h" +#include "MAstroSky2Local.h" #include "MObservatory.h" -#include "MAstroSky2Local.h" -#include "MAstro.h" #include "MLog.h" @@ -94,5 +93,4 @@ return kFALSE; } - fIniSrcPos = *fSrcPos; // look for/create MDCA @@ -181,11 +179,12 @@ // de-rotate the source position - Double_t rotationAngle = fObservatory->RotationAngle(fRA,fDEC,eventTime); + const MAstroSky2Local Observation(eventTime, *fObservatory); + Double_t rotationAngle = Observation.RotationAngle(fRA,fDEC); Float_t c = TMath::Cos(rotationAngle); Float_t s = TMath::Sin(rotationAngle); // perform a rotation of -rotationAngle to move the source back to the "initial" position - Float_t newX = c*fIniSrcPos.GetX()+s*fIniSrcPos.GetY(); - Float_t newY = -s*fIniSrcPos.GetX()+c*fIniSrcPos.GetY(); + Float_t newX = c*fSrcPos->GetX()+s*fSrcPos->GetY(); + Float_t newY = -s*fSrcPos->GetX()+c*fSrcPos->GetY(); #ifdef DEBUG Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.h (revision 3972) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.h (revision 3973) @@ -27,5 +27,4 @@ TString fDCAName; - MSrcPosCam fIniSrcPos; // Copy of the initial source position Double_t fRA; // [rad] Right ascention Double_t fDEC; // [rad] Declination Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcTranslate.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcTranslate.cc (revision 3973) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcTranslate.cc (revision 3973) @@ -0,0 +1,133 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Rico 04/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSrcTranslate +// +// Task to translate the position of the source, starting from previous source +// position (fTranslationIsRelative=kTRUE) or from camera center +// (fTranslationIsRelative=kFALSE) +// +// Input Containers: +// MSrcPosCam +// +// Output Containers: +// MSrcPosCam +// MDCA +// +////////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include "MParList.h" + +#include "MSrcTranslate.h" +#include "MSrcPosCam.h" +#include "MDCA.h" + + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSrcTranslate); + +using namespace std; + +static const TString gsDefName = "MSrcTranslate"; +static const TString gsDefTitle = "De-rotate position of the source"; + +// ------------------------------------------------------------------------- +// +// Default constructor. The first (second) argument is the name of a container +// containing the source position in the camera plain, MScrPosCam (MDCA). +// The default is "MSrcPosCam" ("MDCA"). +// +MSrcTranslate::MSrcTranslate(const char* srcPos, const char* dca, const char *name, const char *title) + : fSrcPos(NULL), fDCA(NULL), fShiftX(0.), fShiftY(0.), fTranslationIsRelative(kTRUE) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + fSrcPosName = srcPos; + fDCAName = dca; +} + +// ------------------------------------------------------------------------- +// +// Look for/create the needed containers +// +Int_t MSrcTranslate::PreProcess(MParList *pList) +{ + // look for/create MSrcPosCam + fSrcPos = (MSrcPosCam*)pList->FindObject(AddSerialNumber(fSrcPosName), "MSrcPosCam"); + if (!fSrcPos) + { + *fLog << warn << AddSerialNumber(fSrcPosName) << " [MSrcPosCam] not found... creating default container." << endl; + fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam", AddSerialNumber(fSrcPosName)); + if(!fSrcPos) + return kFALSE; + } + + // look for/create MDCA + fDCA = (MDCA*)pList->FindObject(AddSerialNumber(fDCAName), "MDCA"); + if (!fDCA) + { + *fLog << warn << AddSerialNumber(fDCAName) << " [MDCA] not found... creating default container." << endl; + fDCA = (MDCA*)pList->FindCreateObj("MDCA", AddSerialNumber(fDCAName)); + if(!fDCA) + return kFALSE; + } + + if(fShiftX==0. && fShiftY==0.) + *fLog << warn << "MSrcTranslate::PreProcess Warning: Null translation" << endl; + + return kTRUE; +} + +// ------------------------------------------------------------------------- +// +// Perform the translation of the source position +// +Int_t MSrcTranslate::Process() +{ + Double_t newX,newY; + + if(fTranslationIsRelative) + { + newX=fSrcPos->GetX()+fShiftX; + newY=fSrcPos->GetY()+fShiftY; + } + else + { + newX=fShiftX; + newY=fShiftY; + } + + fSrcPos->SetX(newX); + fSrcPos->SetY(newY); + fDCA->SetRefPoint(newX,newY); + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcTranslate.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcTranslate.h (revision 3973) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcTranslate.h (revision 3973) @@ -0,0 +1,44 @@ +#ifndef MARS_MSrcTranslate +#define MARS_MSrcTranslate + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#include +#include "MTime.h" +#include "MSrcPosCam.h" + +class MDCA; +class MObservatory; +class MRawEvtHeader; +class MRawRunHeader; + +class MSrcTranslate : public MTask +{ +private: + MSrcPosCam* fSrcPos; //! Pointer to the source position + MDCA* fDCA; //! Pointer to the MDCA object + + TString fSrcPosName; + TString fDCAName; + + Double_t fShiftX; + Double_t fShiftY; + Bool_t fTranslationIsRelative; + + Int_t PreProcess(MParList *plist); + Int_t Process(); + +public: + MSrcTranslate(const char* src="MSrcPosCam", const char* dca="MDCA", + const char* name=NULL, const char* title=NULL); + + void SetTranslation(Double_t x=0.,Double_t y=0) {fShiftX=x,fShiftY=y;}; + void SetRelativeTranslation(Bool_t inp=kTRUE) {fTranslationIsRelative=inp;}; + + ClassDef(MSrcTranslate, 0) // task to rotate the position of the source as a function of Azimuth and Zenith angles +}; + +#endif + Index: trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile (revision 3972) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile (revision 3973) @@ -49,6 +49,7 @@ MPSFFit.cc \ MPSFFitCalc.cc \ - MSrcPosFromFile.cc -# MSrcRotate.cc + MSrcPosFromFile.cc \ + MSrcRotate.cc \ + MSrcTranslate.cc Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/signal.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/signal.C (revision 3973) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/signal.C (revision 3973) @@ -0,0 +1,142 @@ +///////////////////////////////////////////////////////////////////////////// +// +// This macro makes the two-dimension plot of the False Source Method: +// root > .x signal.C(TString name) +// +// Ester Aliu +// Oscar Blanch +// Javier Rico +//////////////////////////////////////////////////////////////////////////// +void signal(TString fname ="hillasCrab/falseSourceCrab20040215Derotate.root") +{ + Float_t excess; + Float_t significance; + Float_t array[2]; + + Int_t nbin2d = 10; + Int_t ntotal = 2*nbin2d+1; + + Int_t nx = 0; + Int_t ny = 0; + + TH2F *hsignif = new TH2F("hsignif","Significance", ntotal, -1.0, 1.0, ntotal, -1.0, 1.0); + TH2F *hexcess = new TH2F("hexcess","Excess", ntotal, -1.0, 1.0, ntotal, -1.0, 1.0); + + // loop on the histos and compute excess and significance + for( Int_t i = -nbin2d; i <= nbin2d ; i++) + { + nx++; + ny = 0; + for( Int_t j = -nbin2d; j <= nbin2d ; j++) + { + ny++; + step_signal(i,j,&array[0], fname); + + excess = array[0]; + significance = array[1]; + hexcess->SetBinContent(nx, ny, excess); + hsignif->SetBinContent(nx, ny, significance); + } + } + + + // Plot + TPostScript myps("signal.ps",111); + myps.Range(15,15); + TCanvas *c2d = new TCanvas("c2d", "Matrices", 0, 0, 800,800); + c2d->Divide(2,2); + gStyle->SetPalette(1); + gStyle->SetOptStat(kFALSE); + + c2d->cd(1); + hexcess->Draw("colz"); + hexcess->SetXTitle("X position (deg)"); + hexcess->SetYTitle("Y position (deg)"); + gPad->Update(); + + c2d->cd(2); + hsignif->Draw("colz"); + hsignif->SetXTitle("X position (deg)"); + hsignif->SetYTitle("Y position (deg)"); + gPad->Update(); + + c2d->cd(3); + hexcess->Draw("lego2"); + hexcess->SetXTitle("X position (deg)"); + hexcess->SetYTitle("Y position (deg)"); + gPad->Update(); + + c2d->cd(4); + hsignif->Draw("lego2"); + hsignif->SetXTitle("X position (deg)"); + hsignif->SetYTitle("Y position (deg)"); + gPad->Update(); + + myps.Close(); +} + + +void step_signal(Int_t i, Int_t j, Float_t *array, TString fname) +{ + cout << "Bin (" << i << "," << j << ")" << endl; + TH1F *onhisto; + TH1F *offhisto; + + // open file and access to the histograms + TFile *f = new TFile(fname); + Char_t name[20]; + Char_t title[50]; + + // histo name + sprintf(name,"hOnAlpha[%d][%d]", i, j); + sprintf(title,"Alpha-Plot(On data) (%d ,%d)", i, j); + onhisto = (TH1F*)gDirectory->Get(name); + + sprintf(name,"hOffAlpha[%d][%d]", i, j); + sprintf(title,"Alpha-Plot(Off data) (%d ,%d)", i, j); + offhisto = (TH1F*)gDirectory->Get(name); + + // line/normalization + const Int_t nbins = 18; + const Int_t inibin = 20./90.*nbins+1; + Float_t level=0; + Float_t leveloff=0; + Float_t levelofferror=0; + for(Int_t ibin = inibin; ibin<=nbins;ibin++) + { + level+=onhisto->GetBinContent(ibin); + leveloff+=offhisto->GetBinContent(ibin); + } + level/=(nbins-inibin+1); + leveloff/=(nbins-inibin+1); + + // normalize on/off + offhisto->Sumw2(); // needed to compute correct errors after normalization + const Float_t norm = level/leveloff; + cout << "Normalizing by factor " << norm <Scale(norm); + + // significance: + Float_t sig=0,bg=0,esig=0,ebg=0; + Float_t significance=0; + const Int_t signbins = inibin-1; + for(Int_t ibin = 1; ibin<=signbins;ibin++) + { + sig += onhisto->GetBinContent(ibin); + esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin); + bg += offhisto->GetBinContent(ibin); + ebg += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin); + } + Float_t error= TMath::Sqrt(esig+ebg); + Float_t excess = sig-bg; + Float_t significance = excess/error; + + cout << "Excess: " << excess << endl; + cout << "Significance: " << significance << endl; + + array[0] = excess; + array[1] = significance; +} + + + Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/signalPoint.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/signalPoint.C (revision 3973) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/signalPoint.C (revision 3973) @@ -0,0 +1,131 @@ +///////////////////////////////////////////////////////////////////////////// +// +// This macro makes the alpha-plot for one of the positions of the +// False Source Method: +// +// Ester Aliu +// Oscar Blanch +// Javier Rico +//////////////////////////////////////////////////////////////////////////// + +void signalPoint(Int_t i=0, Int_t j=0, TString fname="hillasCrab/falseSourceCrab20040215.root") +{ + // constants + const Int_t nbins = 18; + + // general settings + gROOT->Reset(); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetFrameBorderMode(0); + gStyle->SetOptTitle(0); + gStyle->SetTitleOffset(1.7,"y"); + gStyle->SetPadLeftMargin(0.15); + gStyle->SetOptStat(111110); + gStyle->SetStatColor(0); + gStyle->SetStatBorderSize(1); + gStyle->SetStatW(0.2); + gStyle->SetStatH(0.1); + gStyle->SetStatX(0.9); + gStyle->SetStatY(0.9); + + + // open file and access to the histograms + TFile *f = new TFile(fname); + + Char_t name[20]; + Char_t title[50]; + + sprintf(name,"hOnAlpha[%d][%d]", i, j); + sprintf(title,"Alpha-Plot(On data) (%d ,%d)", i, j); + TH1F* onhisto = (TH1F*)gDirectory->Get(name); + + sprintf(name,"hOffAlpha[%d][%d]", i, j); + sprintf(title,"Alpha-Plot(Off data) (%d ,%d)", i, j); + TH1F* offhisto = (TH1F*)gDirectory->Get(name); + + TPostScript myps("signalpoint.ps",111); + myps.Range(15,15); + TCanvas *c1 = new TCanvas("c1","c1",800,600); + c1->cd(1); + + // line/normalization + const Int_t inibin = 20./90.*nbins+1; + Float_t level=0; + Float_t leveloff=0; + Float_t levelofferror=0; + for(Int_t ibin = inibin; ibin<=nbins;ibin++) + { + level+=onhisto->GetBinContent(ibin); + leveloff+=offhisto->GetBinContent(ibin); + } + level/=(nbins-inibin+1); + leveloff/=(nbins-inibin+1); + + // normalize on/off + offhisto->Sumw2(); // needed to compute correct errors after normalization + const Float_t norm = level/leveloff; + cout << "Normalizing by factor " << norm <Scale(norm); + + // significance: + Float_t sig=0,bg=0,esig=0,ebg=0; + Float_t significance=0; + const Int_t signbins = inibin-1; + for(Int_t ibin = 1; ibin<=signbins;ibin++) + { + // Float_t sigma = (onhisto->GetBinContent(ibin)-level)/onhisto->GetBinError(ibin); + // significance+=sigma*sigma; + sig += onhisto->GetBinContent(ibin); + esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin); + bg += offhisto->GetBinContent(ibin); + ebg += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin); + } + Float_t error= TMath::Sqrt(esig+ebg); + significance = (sig-bg)/error; + + cout << "Significance: " << significance << endl; + + // plot + onhisto->SetXTitle("alpha (deg)"); + onhisto->SetYTitle("Entries"); + onhisto->SetMarkerColor(4); + onhisto->SetMarkerStyle(20); + onhisto->SetMarkerSize(.7); + onhisto->SetLineColor(4); + onhisto->SetMaximum(onhisto->GetBinContent(nbins)*1.5); + + offhisto->SetFillColor(2); + offhisto->SetMaximum(offhisto->GetBinContent(nbins)*1.5); + offhisto->SetXTitle("alpha (deg)"); + offhisto->SetYTitle("Entries"); + + offhisto->Draw("HIST"); + offhisto->Draw("ESAME"); + onhisto->Draw("ESAMES"); + // onhisto->Draw("E"); + + // move stat box to make them all visible + gPad->Update(); + TPaveStats* offpavstat = (TPaveStats*) offhisto->GetListOfFunctions()->FindObject("stats"); + if(offpavstat) + { + Float_t shiftx = offpavstat->GetX2NDC()-offpavstat->GetX1NDC(); + offpavstat->SetX1NDC(offpavstat->GetX1NDC()-shiftx); + offpavstat->SetX2NDC(offpavstat->GetX2NDC()-shiftx); + } + + // draw line + TLine* line = new TLine(0.,level,90.,level); + line->SetLineColor(2); + line->SetLineStyle(2); + // line->Draw(); + + gPad->Modified(); + gPad->Update(); + + myps.Close(); +} + + Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/programs/Makefile (revision 3972) +++ trunk/MagicSoft/Mars/mtemp/mifae/programs/Makefile (revision 3973) @@ -18,5 +18,5 @@ include ../../../Makefile.conf.general -PROGRAMS = makeHillas psffit +PROGRAMS = makeHillas psffit falseSource srcPos SOLIB = ../../../libmars.so Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/falseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/programs/falseSource.cc (revision 3973) +++ trunk/MagicSoft/Mars/mtemp/mifae/programs/falseSource.cc (revision 3973) @@ -0,0 +1,476 @@ +//////////////////////////////////////////////////////////////////////////////////// +// +// _____False Source Method macro_____ +// +// Take as input root files with hillas parameters, perform derotation (optional) +// and generates the alpha histograms for On and Off for different positions of the +// source. These histos can be ploted with signal.C and signalPoint.C +// +// Ester Aliu +// Oscar Blanch +// Javier Rico +//////////////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include "TString.h" +#include "TH1F.h" +#include "TFile.h" + +#include "MHillasSrcCalc.h" +#include "MSrcPosCam.h" +#include "MHillasSrc.h" +#include "MSrcRotate.h" +#include "MParList.h" +#include "MTaskList.h" +#include "MHillas.h" +#include "MReadTree.h" +#include "MEvtLoop.h" +#include "MLog.h" +#include "MArgs.h" + +using namespace std; + +Bool_t readDatacards(TString& filename); +void falseSource(); + +//----------------------------------------------------------------------------- +// declaration of variables read from datacards +//----------------------------------------------------------------------------- + +TString onFile; +TString offFile; +TString outputFile; +Bool_t kRotate=1; +Double_t fRA= -1.; +Double_t fDEC= -1.; +ULong_t nmaxevents=999999999; + +// cuts (default values to be overriden by datacards) +Float_t minsize = 0; // minimum size (# of photons) +Float_t maxsize = 9999999; // maximum size (# of photons) +Float_t mindist = 0; // minimum distance cut (deg) +Float_t maxdist = 10; // maximum distance cut (deg) +Float_t minwidth = 0.; // minimum width cut (deg) +Float_t maxwidth = 10; // maximum width cut (deg) +Float_t minlength = 0; // minimum length cut (deg) +Float_t maxlength = 10; // maximum length cut (deg) +Float_t maxXY = maxdist; // maximum X and Y distance from the camera center + +// binning +Int_t binning = 21; // number of bins in false source search (one coordinate) +Int_t nbin = 18; // number of bins for alpha plots +Float_t field = 2.; // width of the examined field (degrees) + +//----------------------------------------------------------------------------- +// constants +//----------------------------------------------------------------------------- + +const TString defaultcard="falsesource.datacard"; + +const Float_t alphamin = 0.; // minimum value in alpha (degrees) +const Float_t alphamax = 90.; // maximum value in alpha (degrees) +const Float_t conver = 189./0.6; // conversion factor degrees to mm + +//----------------------------------------------------------------------------- + +static void Usage() +{ + gLog <" << endl << endl; + gLog << " : datacards file name (dafault falsesource.datacards)" << endl; + gLog << " -?/-h: This help" << endl << endl; +} + +//----------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + // evaluate arguments + MArgs arg(argc, argv); + if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) + { + Usage(); + return -1; + } + + TString datacard = arg.GetArgumentStr(0); + if(!datacard.Length()) + datacard = defaultcard; + + if(!readDatacards(datacard)) + { + cout << "Error reading datacards. Stoping" << endl; + return -1; + } + falseSource(); +} + +void falseSource() +{ + // CUTS + + // variable declaration + Float_t xpos; + Float_t ypos; + Char_t name[20]; + Char_t title[50]; + Float_t stepDegrees = field/(binning-1); // degrees + Float_t step = stepDegrees*conver; // mm + Int_t ival = (Int_t) ceil((Float_t)(binning-1)/2); // index bound + + // Hillas variables + Float_t alpha; + Float_t size; + Float_t dist; + Float_t length; + Float_t width; + Float_t meanX; + Float_t meanY; + + // create the histograms (empty) + TH1F *hOnAlpha[binning][binning]; + TH1F *hOffAlpha[binning][binning]; + for(Int_t i = -ival; i <= ival ; i++){ + for(Int_t j = -ival; j <= ival ; j++){ + + sprintf(name,"hOnAlpha[%d][%d]", i, j); + sprintf(title,"Alpha-Plot(On data) (%d ,%d)", i, j); + hOnAlpha[i+ival][j+ival] = new TH1F(name, title, nbin, alphamin, alphamax); + + sprintf(name,"hOffAlpha[%d][%d]", i, j); + sprintf(title,"Alpha-Plot(Off data) (%d ,%d)", i, j); + hOffAlpha[i+ival][j+ival] = new TH1F(name, title, nbin, alphamin, alphamax); + + } + } + + + + /****************************************************************** + FILL ON-DATA HISTOS + *******************************************************************/ + + // source dependent hillas containers/tasks + MHillasSrcCalc* csrc1[binning][binning]; + MSrcPosCam* source[binning][binning]; + MHillasSrc* hillasSrc[binning][binning]; + MSrcRotate* srcrotate[binning][binning]; + + // normal containers/classes + MParList plist; + MTaskList tlist; + MHillas hillas; + + plist.AddToList(&tlist); + plist.AddToList(&hillas); + + MReadTree read("Parameters", onFile); + read.DisableAutoScheme(); + tlist.AddToList(&read); + + Int_t k,l; + char sourceName[100]; + char hillasSrcName[100]; + + // create the tasks/containers for the different source positions + for(Int_t i=-ival;i<=ival;i++) + { + xpos = step*i; + for(Int_t j=-ival;j<=ival;j++) + { + ypos = step*j; + + // name the different containers + if (i<0 && j<0) + { + k = -i; + l = -j; + sprintf(sourceName, "MSrcPosCam_Min%dMin%d", k, l); + sprintf(hillasSrcName, "MHillasSrc_Min%dMin%d", k, l); + } + + if (i<0 && j>=0) + { + k = -i; + sprintf(sourceName, "MSrcPosCam_Min%d%d", k, j); + sprintf(hillasSrcName, "MHillasSrc_Min%d%d", k, j); + } + + if (i>=0 && j<0) + { + l = -j; + sprintf(sourceName, "MSrcPosCam_%dMin%d", i, l); + sprintf(hillasSrcName, "MHillasSrc_%dMin%d", i, l); + } + + if (i>=0 && j>= 0) + { + sprintf(sourceName, "MSrcPosCam_%d%d", i, j); + sprintf(hillasSrcName, "MHillasSrc_%d%d", i, j); + } + + // include containers in parameter list + source[i+ival][j+ival] = new MSrcPosCam(sourceName); + source[i+ival][j+ival]->SetXY(xpos, ypos); + plist.AddToList(source[i+ival][j+ival]); + source[i+ival][j+ival]->Print(); + + hillasSrc[i+ival][j+ival] = new MHillasSrc(hillasSrcName); + plist.AddToList(hillasSrc[i+ival][j+ival]); + + // define the different tasks and include them in the task list + if(kRotate) + { + srcrotate[i+ival][j+ival] = new MSrcRotate(sourceName); + srcrotate[i+ival][j+ival]->SetRAandDEC(fRA,fDEC); + tlist.AddToList(srcrotate[i+ival][j+ival]); + } + + TString HilName = "MHillas"; + csrc1[i+ival][j+ival] = new MHillasSrcCalc(sourceName, hillasSrcName); + csrc1[i+ival][j+ival]->SetInput(HilName); + + tlist.AddToList(csrc1[i+ival][j+ival]); + + } // loop on j (y coordinate) + } // loop on i (x coordinate) + + + // Eventloop + MEvtLoop evtloop; + evtloop.SetParList(&plist); + // MProgressBar bar; + // evtloop.SetProgressBar(&bar); + if (!evtloop.PreProcess()) + return; + + // select the events passing the cuts + UInt_t nread=0; + while (tlist.Process()) + { + for(Int_t i = -ival; i <= ival ; i++) + for(Int_t j = -ival; j <= ival ; j++) + { + width = hillas.GetWidth()/conver; + length = hillas.GetLength()/conver; + size = hillas.GetSize(); + dist = hillasSrc[i+ival][j+ival]->GetDist()/conver; + meanX = hillas.GetMeanX()/conver; + meanY = hillas.GetMeanY()/conver; + + if (widthmindist && distminsize && meanXGetAlpha(); + hOnAlpha[i+ival][j+ival]->Fill(abs(alpha)); + } + } + if(++nread>nmaxevents) break; + } + + // end + evtloop.PostProcess(); + tlist.PrintStatistics(); + + + /****************************************************************** + FILL OFF-DATA HISTOS + *******************************************************************/ + + MReadTree read2("Parameters", offFile); + read2.DisableAutoScheme(); + + tlist.AddToListBefore(&read2, &read, "All"); + tlist.RemoveFromList(&read); + + if (!evtloop.PreProcess()) + return; + + nread=0; + while (tlist.Process()) + { + for(Int_t i = -ival; i <= ival ; i++) + for(Int_t j = -ival; j <= ival ; j++) + { + width = hillas.GetWidth()/conver; + length = hillas.GetLength()/conver; + size = hillas.GetSize(); + dist = hillasSrc[i+ival][j+ival]->GetDist()/conver; + meanX = hillas.GetMeanX()/conver; + meanY = hillas.GetMeanY()/conver; + + if (widthmindist && distminsize && meanXGetAlpha(); + hOffAlpha[i+ival][j+ival]->Fill(abs(alpha)); + } + } + if(++nread>nmaxevents) break; + } + + evtloop.PostProcess(); + + // Save results + TFile *f = new TFile(outputFile, "RECREATE"); + + for(Int_t i = -ival; i <= ival ; i++){ + for(Int_t j = -ival; j <= ival ; j++){ + hOnAlpha[i+ival][j+ival]->Write(); + hOffAlpha[i+ival][j+ival]->Write(); + } + } + + f->Close(); + +} +//----------------------------------------------------------------------- + +Bool_t readDatacards(TString& filename) +{ + ifstream ifun(filename.Data()); + if(!ifun) + { + cout << "File " << filename << " not found" << endl; + return kFALSE; + } + + TString word; + + while(ifun >> word) + { + // skip comments + if(word[0]=='/' && word[1]=='/') + { + while(ifun.get()!='\n'); // skip line + continue; + } + + // number of events + if(strcmp(word.Data(),"NEVENTS")==0) + ifun >> nmaxevents; + + // on-data file name + if(strcmp(word.Data(),"ONFILES")==0) + { + if(onFile.Length()) + cout << "readDataCards Warning: overriding on-data file name" << endl; + ifun >> onFile; + } + + // off-data file name + if(strcmp(word.Data(),"OFFFILES")==0) + { + if(offFile.Length()) + cout << "readDataCards Warning: overriding off-data file name" << endl; + ifun >> offFile; + } + + // output file name + if(strcmp(word.Data(),"OUTFILE")==0) + { + if(outputFile.Length()) + cout << "readDataCards Warning: overriding output file name" << endl; + ifun >> outputFile; + } + + // rotation flag + if(strcmp(word.Data(),"ROTFLAG")==0) + ifun >> kRotate; + + // source celestial coordinates + if(strcmp(word.Data(),"SRCCOORDS")==0) + { + ifun >> fRA; + ifun >> fDEC; + } + + // field width + if(strcmp(word.Data(),"FIELD")==0) + ifun >> field; + + // binning + if(strcmp(word.Data(),"BINNING")==0) + ifun >> binning; + + + // Number of bins in alpha plots + if(strcmp(word.Data(),"NBIN")==0) + ifun >> nbin; + + + + + // size cut + if(strcmp(word.Data(),"SIZECUT")==0) + { + ifun >> minsize; + ifun >> maxsize; + } + + // dist cut + if(strcmp(word.Data(),"DISTCUT")==0) + { + ifun >> mindist; + ifun >> maxdist; + } + + // width cut + if(strcmp(word.Data(),"WIDTHCUT")==0) + { + ifun >> minwidth; + ifun >> maxwidth; + } + + // length cut + if(strcmp(word.Data(),"LENGTHCUT")==0) + { + ifun >> minlength; + ifun >> maxlength; + } + + // maxX and maxY upper cut + if(strcmp(word.Data(),"CENTERCUT")==0) + ifun >> maxXY; + } + + // check compulsory values + if(!onFile.Length()) + { + cout << "No on-data file name specified" << endl; + return kFALSE; + } + if(!offFile.Length()) + { + cout << "No off-data file name specified" << endl; + return kFALSE; + } + if(!outputFile.Length()) + { + cout << "No output file name specified" << endl; + return kFALSE; + } + + + // Dump read values + cout << "************************************************" << endl; + cout << "* Datacards read from file " << filename << endl; + cout << "************************************************" << endl; + cout << "Maximum number of input events: " << nmaxevents << endl; + cout << "On-data file name(s): " << onFile << endl; + cout << "Off-data file name(s): " << offFile << endl; + cout << "Output file name: " << outputFile << endl; + cout << "De-rotation flag " << kRotate << endl; + cout << "Source celestial coordiantes (rad): RA = " << fRA << ", DEC = " << fDEC << endl; + cout << "Field width (degrees): " << field << endl; + cout << "Field binning: " << binning << endl; + cout << "Number of alpha plot bins: " << nbin << endl; + cout << "Size cuts (# of photons): ("< +// Oscar Blanch +// Javier Rico +//////////////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include "TString.h" + +#include "MHillasSrcCalc.h" +#include "MHillasSrc.h" +#include "MSrcRotate.h" +#include "MSrcTranslate.h" +#include "MParList.h" +#include "MTaskList.h" +#include "MHillas.h" +#include "MReadTree.h" +#include "MEvtLoop.h" +#include "MLog.h" +#include "MArgs.h" +#include "MWriteRootFile.h" + +using namespace std; + +Bool_t readDatacards(TString& filename); +void srcPos(); + +//----------------------------------------------------------------------------- +// declaration of variables read from datacards +//----------------------------------------------------------------------------- + +TString inputFile; +TString outputFile; +ULong_t nmaxevents=999999999; +Float_t xsrcpos=0.; +Float_t ysrcpos=0.; +Bool_t kRotate=1; +Bool_t kSrcPolicy=kTRUE; +Double_t fRA= -1.; +Double_t fDEC= -1.; + +//----------------------------------------------------------------------------- +// constants +//----------------------------------------------------------------------------- + +const TString defaultcard="srcpos.datacard"; + +const Float_t alphamin = 0.; // minimum value in alpha (degrees) +const Float_t alphamax = 90.; // maximum value in alpha (degrees) +const Float_t conver = 189./0.6; // conversion factor degrees to mm + +//----------------------------------------------------------------------------- + +static void Usage() +{ + gLog <" << endl << endl; + gLog << " : datacards file name (dafault " << defaultcard <<")" << endl; + gLog << " -?/-h: This help" << endl << endl; +} + +//----------------------------------------------------------------------------- +int main(int argc, char **argv) +{ + // evaluate arguments + MArgs arg(argc, argv); + if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) + { + Usage(); + return -1; + } + + TString datacard = arg.GetArgumentStr(0); + if(!datacard.Length()) + datacard = defaultcard; + + if(!readDatacards(datacard)) + { + cout << "Error reading datacards. Stoping" << endl; + return -1; + } + srcPos(); +} + +//----------------------------------------------------------------------------- +void srcPos() +{ + // variable declaration + Float_t xpos=xsrcpos*conver; // [mm] + Float_t ypos=ysrcpos*conver; // [mm] + + // containers + MParList plist; + MTaskList tlist; + + // include containers in parameter list + plist.AddToList(&tlist); + + // tasks + MReadTree read("Parameters", inputFile); + read.DisableAutoScheme(); + + MSrcTranslate srctranslate; + srctranslate.SetTranslation(xpos,ypos); + srctranslate.SetRelativeTranslation(kSrcPolicy); + + MSrcRotate srcrotate; + srcrotate.SetRAandDEC(fRA,fDEC); + + MHillasSrcCalc csrc1; + + MWriteRootFile write(outputFile,"RECREATE"); + write.AddContainer("MHillas" , "Parameters"); + write.AddContainer("MHillasSrc" , "Parameters"); + write.AddContainer("MHillasExt" , "Parameters"); + write.AddContainer("MNewImagePar" , "Parameters"); + write.AddContainer("MRawEvtHeader" , "Parameters"); + write.AddContainer("MRawRunHeader" , "Parameters"); + write.AddContainer("MConcentration" , "Parameters"); + write.AddContainer("MSrcPosCam" , "Parameters"); + + // include tasks in task list + tlist.AddToList(&read); + tlist.AddToList(&srctranslate); + if(kRotate) + tlist.AddToList(&srcrotate); + tlist.AddToList(&csrc1); + tlist.AddToList(&write); + + // Eventloop + MEvtLoop evtloop; + evtloop.SetParList(&plist); + if (!evtloop.Eventloop(nmaxevents)) + return; + + tlist.PrintStatistics(); + +} +//----------------------------------------------------------------------- + +Bool_t readDatacards(TString& filename) +{ + ifstream ifun(filename.Data()); + if(!ifun) + { + cout << "File " << filename << " not found" << endl; + return kFALSE; + } + + TString word; + + while(ifun >> word) + { + // skip comments + if(word[0]=='/' && word[1]=='/') + { + while(ifun.get()!='\n'); // skip line + continue; + } + + // number of events + if(strcmp(word.Data(),"NEVENTS")==0) + ifun >> nmaxevents; + + // input file name + if(strcmp(word.Data(),"INPUTFILES")==0) + { + if(inputFile.Length()) + cout << "readDataCards Warning: overriding on-data file name" << endl; + ifun >> inputFile; + } + + // output file name + if(strcmp(word.Data(),"OUTFILE")==0) + { + if(outputFile.Length()) + cout << "readDataCards Warning: overriding output file name" << endl; + ifun >> outputFile; + } + + // source position + if(strcmp(word.Data(),"SRCPOS")==0) + { + ifun >> xsrcpos; + ifun >> ysrcpos; + } + + // source celestial coordinates + if(strcmp(word.Data(),"SRCCOORDS")==0) + { + ifun >> fRA; + ifun >> fDEC; + } + + // source movement policy flag + if(strcmp(word.Data(),"SRCABS")==0) + ifun >> kSrcPolicy; + + // rotation flag + if(strcmp(word.Data(),"ROTFLAG")==0) + ifun >> kRotate; + } + + // check compulsory values + if(!inputFile.Length()) + { + cout << "No on-data file name specified" << endl; + return kFALSE; + } + + if(!outputFile.Length()) + { + cout << "No output file name specified" << endl; + return kFALSE; + } + + if(xsrcpos==0 && ysrcpos==0) + { + cout << "Source position is center of the camera (as in input file)" << endl; + return kFALSE; + } + + // Dump read values + cout << "************************************************" << endl; + cout << "* Datacards read from file " << filename << endl; + cout << "************************************************" << endl; + cout << "Maximum number of input events: " << nmaxevents << endl; + cout << "Input file name(s): " << inputFile << endl; + cout << "Output file name: " << outputFile << endl; + cout << "Source position (degrees) X=" << xsrcpos << ", Y="<