Index: /trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 4049) +++ /trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 4050) @@ -18,6 +18,16 @@ -*-*- END OF LINE -*-*- + 2004/05/11: Javier Rico + * programs/optimizeCuts.cc, programs/optimizecuts.datacard + - added -2004/05/11: Ester Aliu + * programs/Makefile + - include optimizeCuts + + * programs/srcPos.cc, programs/srcpos.datacard + - add reference rotation angle + + + 2004/05/11: Ester Aliu * macros/signal.cc - change some things @@ -27,9 +37,9 @@ - add library -2004/05/11: Ester Aliu + 2004/05/11: Ester Aliu * programs/makeHillas.cc, programs/makehillas.datacard - add island calculations -2004/05/11: Javier Rico + 2004/05/11: Javier Rico * programs/makeHillas.cc - prevent it from bombing because of the mjobs Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc (revision 4049) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc (revision 4050) @@ -180,5 +180,6 @@ // de-rotate the source position const MAstroSky2Local Observation(eventTime, *fObservatory); - Double_t rotationAngle = Observation.RotationAngle(fRA,fDEC); + const MAstroSky2Local RefObservation(fRefMJD, *fObservatory); + Double_t rotationAngle = Observation.RotationAngle(fRA,fDEC)-RefObservation.RotationAngle(fRA,fDEC); Float_t c = TMath::Cos(rotationAngle); Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.h (revision 4049) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.h (revision 4050) @@ -27,12 +27,13 @@ TString fDCAName; - Double_t fRA; // [rad] Right ascention - Double_t fDEC; // [rad] Declination - UInt_t fNEvts; // Number of events in file - UInt_t fFirstEvt; // Event index for first event in file - MTime fIniTime; // Run initial time - MTime fFinTime; // Run final time - Double_t fDeltaT; // DeltaT between two events - UInt_t fRunNumber; // Current run number + Double_t fRA; // [rad] Right ascention + Double_t fDEC; // [rad] Declination + Double_t fRefMJD; // [MJ date] reference time for rotation + UInt_t fNEvts; // Number of events in file + UInt_t fFirstEvt; // Event index for first event in file + MTime fIniTime; // Run initial time + MTime fFinTime; // Run final time + Double_t fDeltaT; // DeltaT between two events + UInt_t fRunNumber; // Current run number Int_t PreProcess(MParList *plist); @@ -43,5 +44,5 @@ const char* name=NULL, const char* title=NULL); - void SetRAandDEC(Double_t ra, Double_t dec) {fRA=ra;fDEC=dec;}; + void SetRAandDECandRefMJD(Double_t ra, Double_t dec, Double_t ref) {fRA=ra;fDEC=dec;fRefMJD=ref;} ClassDef(MSrcRotate, 0) // task to rotate the position of the source as a function of Azimuth and Zenith angles Index: /trunk/MagicSoft/Mars/mtemp/mifae/programs/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/programs/Makefile (revision 4049) +++ /trunk/MagicSoft/Mars/mtemp/mifae/programs/Makefile (revision 4050) @@ -18,5 +18,5 @@ include ../../../Makefile.conf.general -PROGRAMS = makeHillas psffit falseSource srcPos +PROGRAMS = makeHillas psffit falseSource srcPos optimizeCuts SOLIB = $(MARSSYS)/libmars.so Index: /trunk/MagicSoft/Mars/mtemp/mifae/programs/makehillas.datacard =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/programs/makehillas.datacard (revision 4049) +++ /trunk/MagicSoft/Mars/mtemp/mifae/programs/makehillas.datacard (revision 4050) @@ -1,5 +1,5 @@ // Maximun number of (data) events to be processed) -NEVENTS 9999999 +NEVENTS 10 // data file directory Index: /trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizeCuts.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizeCuts.cc (revision 4050) +++ /trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizeCuts.cc (revision 4050) @@ -0,0 +1,326 @@ +//////////////////////////////////////////////////////////////////////////////////// +// +// _____ Optimize cuts _____ +// +// Take as input root files with hillas parameters and compute the optimal cuts by +// mapping the width/length plane (only upper cuts so far) +// +// Jose Flix +// Javier Rico +//////////////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include "TString.h" +#include "TChain.h" +#include "TH1F.h" + +#include "MArgs.h" +#include "MLog.h" + +using namespace std; + +Bool_t readDatacards(TString& filename); +void optimizeCuts(); + + +//----------------------------------------------------------------------------- +// declaration of variables read from datacards +//----------------------------------------------------------------------------- + +TString onFile; +TString offFile; +TString outputFile; +ULong_t nmaxevents=999999999; +Float_t lowdistcut=0; +Float_t uppdistcut=10; +Float_t sizecut=1800; +UInt_t onRate=1; +Float_t lowlgcut=0.01; +Float_t upplgcut=0.3; +Float_t lgstep=0.005; +Float_t lowwdcut=0.01; +Float_t uppwdcut=0.3; +Float_t wdstep=0.005; + +//----------------------------------------------------------------------------- +// constants +//----------------------------------------------------------------------------- + +const TString defaultcard="optimizecuts.datacard"; +const Float_t conver = 189./0.6; // conversion factor degrees to mm +const Int_t nbins=18; //number of bins in alpha plots +const Float_t frontiere=20.; // [deg] value below (above) wich signal (background is computed (normalized) + +//----------------------------------------------------------------------------- + +static void Usage() +{ + gLog <" << endl << endl; + gLog << " : datacards file name (default " << 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; + } + optimizeCuts(); +} + +//----------------------------------------------------------------------------- +void optimizeCuts() +{ + // create the TChains for ON and OFF data + TChain* ton= new TChain("Parameters"); + TChain* toff= new TChain("Parameters"); + if(!ton->Add(onFile)) + { + cout << "On-data file not found or not valid format" << endl; + return; + } + if(!toff->Add(offFile)) + { + cout << "Off-data file not found or not valid format" << endl; + return; + } + + ofstream fout; + fout.open(outputFile.Data()); + + // define aliases + ton->SetAlias("length","MHillas.fLength*0.6/189"); + ton->SetAlias("width","MHillas.fWidth*0.6/189"); + ton->SetAlias("dist","MHillasSrc.fDist*0.6/189"); + ton->SetAlias("conc","MNewImagePar.fConc"); + ton->SetAlias("size","MHillas.fSize"); + ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber"); + ton->SetAlias("alpha","abs(MHillasSrc.fAlpha)"); + + toff->SetAlias("length","MHillas.fLength*0.6/189"); + toff->SetAlias("width","MHillas.fWidth*0.6/189"); + toff->SetAlias("dist","MHillasSrc.fDist*0.6/189"); + toff->SetAlias("conc","MNewImagePar.fConc"); + toff->SetAlias("size","MHillas.fSize"); + toff->SetAlias("event","MRawEvtHeader.fDAQEvtNumber"); + toff->SetAlias("alpha","abs(MHillasSrc.fAlpha)"); + + + // general cut + char gencut[256]; + char evecut[256]; + sprintf(evecut,"event%%%d==0",onRate); + sprintf(gencut,"size>%f && dist>%f && dist<%f",sizecut,lowdistcut,uppdistcut); + cout << "General cut " << gencut << "(" << evecut << ") for on-data" << endl; + + // loop on widht and length cuts + Float_t length=lowlgcut; + while(length<=upplgcut) + { + Float_t width=lowlgcut; + while(width<=uppwdcut) + { + TH1F* onhisto = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90); + TH1F* offhisto = new TH1F("OffHist","Alpha Plot",nbins,0,90); + + // define the width/length cut + char varcut[256]; + sprintf(varcut,"length<%f && width<%f",length,width); + cout << "Cutting " << varcut << endl; + + // define the on/off data cut + char offcut[1024]; + sprintf(offcut,"%s && %s",gencut,varcut); + char oncut[1024]; + sprintf(oncut,"%s && %s",evecut,offcut); + + // Fill the histos + ton->Draw("alpha>>OnHist",oncut); + toff->Draw("alpha>>OffHist",offcut); + + // Normalize histos + const Int_t inibin = (Int_t)(frontiere/90.*nbins+1); + Float_t level=0; + Float_t leveloff=0; + for(Int_t ibin = inibin; ibin<=nbins;ibin++) + { + level+=onhisto->GetBinContent(ibin); + leveloff+=offhisto->GetBinContent(ibin); + } + offhisto->Sumw2(); // needed to compute correct errors after normalization + const Float_t norm = leveloff ? level/leveloff: 1; + cout << "Normalizing by factor " << norm <Scale(norm); + + // compute significance/excess + 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); + significance = error ? (sig-bg)/error : 0; + + cout << "Excess: " << sig-bg << endl; + cout << "N bkg: " << bg << endl; + cout << "Significance: " << significance << endl; + + // save the result in file + fout << width << '\t' + << length << '\t' + << sig-bg << '\t' + << significance << '\t' << endl; + + delete onhisto; + delete offhisto; + + width+=wdstep; + } + length+=lgstep; + } + + fout.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; + + // number of events + if(strcmp(word.Data(),"ONRATE")==0) + ifun >> onRate; + + // input file name (on data) + if(strcmp(word.Data(),"ONFILES")==0) + { + if(onFile.Length()) + cout << "readDataCards Warning: overriding on-data file name" << endl; + ifun >> onFile; + } + + // input file name (off data) + 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; + } + + // size cut + if(strcmp(word.Data(),"SIZECUT")==0) + ifun >> sizecut; + + // dist cuts + if(strcmp(word.Data(),"DISTCUTS")==0) + { + ifun >> lowdistcut ; + ifun >> uppdistcut ; + } + + // length cuts + if(strcmp(word.Data(),"LENGTHCUTS")==0) + { + ifun >> lowlgcut; + ifun >> upplgcut; + ifun >> lgstep; + } + + // width cuts + if(strcmp(word.Data(),"WIDTHCUTS")==0) + { + ifun >> lowwdcut; + ifun >> uppwdcut; + ifun >> wdstep; + } + + } + + // 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 acceptance rate: 1/" << onRate << endl; + cout << "On-data input file name(s): " << onFile << endl; + cout << "Off-data input file name(s): " << offFile << endl; + cout << "Output file name: " << outputFile << endl; + cout << "Size cut: " << sizecut << endl; + cout << "Dist cuts (deg): [" << lowdistcut<< ","<> xsrcpos; ifun >> ysrcpos; + ifun >> mjdpos; } Index: /trunk/MagicSoft/Mars/mtemp/mifae/programs/srcpos.datacard =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/programs/srcpos.datacard (revision 4049) +++ /trunk/MagicSoft/Mars/mtemp/mifae/programs/srcpos.datacard (revision 4050) @@ -4,14 +4,14 @@ // Input file name pattern -INPUTFILES ./srposPrueba.root +INPUTFILES /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OnB.root // output file name -OUTFILE ./srposPrueba2.root +OUTFILE ./srcPosPrueba.root -// X and Y position of the source in degrees -SRCPOS 0.2 -0.1 +// X and Y position of the source in degrees for a MJ date (days) -important only in case of rotation +SRCPOS 0.21 0.176 53050.0773 -// Flag to determine wheter source position is absolute (0) or relative to previous position (1) -SRCABS 1 +// Flag to determine whether source position is absolute (0) or relative to previous position (1) +SRCABS 0 // rotation flag: @@ -21,4 +21,5 @@ // source coordinates (RA DEC in rads) -SRCCOORDS 1.46 0.384 +// SRCCOORDS 1.46 0.384 // (Crab) +SRCCOORDS 2.899 0.667 // (Mrk 421)