Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/programs/Makefile (revision 3965) +++ 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="<