Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 4056) @@ -18,4 +18,19 @@ -*-*- END OF LINE -*-*- + + + + 2004/05/12: Javier Rico + * macros/alpha.C, macros/plotOptimal.C + - added + + * library/MSrcRotate.[cc.h] + - correct bug in computation of run time for rotation + + * programs/srcPos.cc + - change mjdpos type to Double_t + + + 2004/05/11: Javier Rico * programs/optimizeCuts.cc, programs/optimizecuts.datacard Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc (revision 4056) @@ -149,5 +149,4 @@ // save the number of events, initial and final times fNEvts = fRunHeader->GetNumEvents(); - fFirstEvt = fEvtHeader->GetDAQEvtNumber(); fIniTime = fRunHeader->GetRunStart(); fFinTime = fRunHeader->GetRunEnd(); @@ -164,5 +163,4 @@ cout << "Run number: " << fRunHeader->GetRunNumber() << endl; cout << "Number of events: " << fNEvts << endl; - cout << "First event: " << fFirstEvt << endl; cout << "Initial MJD date: " << fIniTime.GetMjd() << endl; cout << "Final MJD date: " << fFinTime.GetMjd() << endl; @@ -175,10 +173,22 @@ // FIXME: for the time being, this is computed by assuming constant event rate MTime eventTime; - Double_t newMJD = fIniTime.GetMjd() + (fFinTime.GetMjd()-fIniTime.GetMjd())*(fEvtHeader->GetDAQEvtNumber()-fFirstEvt)/fNEvts; - eventTime.SetMjd(newMJD); - + Double_t newMJD = fIniTime.GetMjd() + (fFinTime.GetMjd()-fIniTime.GetMjd())*fEvtHeader->GetDAQEvtNumber()/fNEvts; + eventTime.SetMjd(newMJD); + MTime refTime; + refTime.SetMjd(fRefMJD); + // de-rotate the source position const MAstroSky2Local Observation(eventTime, *fObservatory); - const MAstroSky2Local RefObservation(fRefMJD, *fObservatory); + const MAstroSky2Local RefObservation(refTime, *fObservatory); + +#ifdef DEBUG + printf("Run:%d, Event:%d, iniMJD=%15.5f, finMJD=%15.5f, fDeltaT=%15.5f, newMJD=%15.5f, fRefMJD=%15.5f, rotation=%15.5f, ref=%15.5f\n", + fRunHeader->GetRunNumber(),fEvtHeader->GetDAQEvtNumber(), + fIniTime.GetMjd(),fFinTime.GetMjd(),fDeltaT, + newMJD,fRefMJD,Observation.RotationAngle(fRA,fDEC), + RefObservation.RotationAngle(fRA,fDEC)); + // cout << "newMJD=" << newMJD << ", fRefMJD="<GetX()+s*fSrcPos->GetY(); - Float_t newY = -s*fSrcPos->GetX()+c*fSrcPos->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 4055) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.h (revision 4056) @@ -18,9 +18,9 @@ { private: - MSrcPosCam* fSrcPos; //! Pointer to the source position - MDCA* fDCA; //! Pointer to the MDCA object - MRawEvtHeader* fEvtHeader; //! Pointer to the event header - MRawRunHeader* fRunHeader; //! Pointer to the run header - MObservatory* fObservatory; //! Pointer to the MObservatory + MSrcPosCam* fSrcPos; // Pointer to the source position + MDCA* fDCA; // Pointer to the MDCA object + MRawEvtHeader* fEvtHeader; // Pointer to the event header + MRawRunHeader* fRunHeader; // Pointer to the run header + MObservatory* fObservatory; // Pointer to the MObservatory TString fSrcPosName; @@ -31,5 +31,4 @@ 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 Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/alpha.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/alpha.C (revision 4056) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/alpha.C (revision 4056) @@ -0,0 +1,211 @@ +/* ======================================================================== *\! +! Author(s): Oscar Blanch & Javier Rico +\* ======================================================================== */ + +void alpha( + Float_t cutsize=2000, + Float_t upcutsize=9999999, + Float_t lengthlowcut=0., + Float_t lengthcut=0.26, + Float_t widthlowcut=0., + Float_t widthcut=0.12 +) +{ + // constants + const Int_t nbins = 18; + TString psname("alphaMrk.gif"); + // 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); + + // define files, chain, trees, canvas... + TChain* ton= new TChain("Parameters"); + + // ON data + // ton->Add("/disc01/HillasFiles/crabOnMisp3015_A.root"); + // ton->Add("/disc01/HillasFiles/crabOnMisp3015_B.root"); + // ton->Add("/disc01/HillasFiles/crabOnMisp3015_C.root"); + // ton->Add("/disc01/HillasFiles/crabOnMisp3015_D.root"); + // ton->Add("hillasCrab/crab20040215OnA.root"); + // ton->Add("hillasCrab/crab20040215OnB.root"); + // ton->Add("hillasCrab/crab20040215OnC.root"); + // ton->Add("hillasCrab/crab20040215OnD.root"); + // ton->Add("hillasCrab/crab20040215OnRotate2_-1.root"); + ton->Add("/disc01/HillasFiles/Mrk421/mrk421OnMisp3015_A.root"); + ton->Add("/disc01/HillasFiles/Mrk421/mrk421OnMisp3015_B.root"); + //ton->Add("/disc01/HillasFiles/Mrk421/mrk421OnMisp3015_C.root"); + // ton->Add("~aliu/MagicData/CrabNebula/2004_03_19/crabOn3015_A_-1-2.root"); + // ton->Add("~aliu/MagicData/CrabNebula/2004_03_19/crabOn3015_B_-1-1.root"); + //ton->Add("~aliu/MagicData/CrabNebula/2004_03_19/crabOn3015_C.root"); + + TChain* toff= new TChain("Parameters"); + + // OFF data + //toff->Add("/disc01/HillasFiles/crabOffMisp3015_A.root"); + // toff->Add("/disc01/HillasFiles/crabOffMisp3015_B.root"); + // toff->Add("/disc01/HillasFiles/crabOffMisp3015_C.root"); + // toff->Add("/disc01/HillasFiles/crabOffMisp3015_D.root"); + // toff->Add("/disc01/HillasFiles/crabOffMisp3015_E.root"); + // toff->Add("/disc01/HillasFiles/crabOffMisp3015_F.root"); + // toff->Add("/disc01/HillasFiles/crabOffMisp3015_G.root"); + // toff->Add("/disc01/HillasFiles/crabOffMisp3015_H.root"); + // toff->Add("hillasCrab/crab20040215OffA.root"); + // toff->Add("hillasCrab/crab20040215OffB.root"); + // toff->Add("hillasCrab/crab20040215OffC.root"); + // toff->Add("hillasCrab/crab20040215OffD.root"); + // toff->Add("hillasCrab/crab20040215OffE.root"); + // toff->Add("hillasCrab/crab20040215OffF.root"); + // toff->Add("hillasCrab/crab20040215OffG.root"); + // toff->Add("hillasCrab/crab20040215OffH.root"); + // toff->Add("hillasCrab/crab20040215OffRotate2_-1.root"); + toff->Add("/disc01/HillasFiles/OffMrk421/mrk421OffMisp3015_A.root"); + toff->Add("/disc01/HillasFiles/OffMrk421/mrk421OffMisp3015_B.root"); + toff->Add("/disc01/HillasFiles/OffMrk421/mrk421OffMisp3015_C.root"); + // toff->Add("~aliu/MagicData/CrabNebula/2004_03_19/crabOff3015_-1-1.root"); + // toff->Add("~aliu/MagicData/CrabNebula/2004_03_19/crabOff3015_A.root"); + // toff->Add("~aliu/MagicData/CrabNebula/2004_03_19/crabOff3015_B.root"); + // toff->Add("~aliu/MagicData/CrabNebula/2004_03_19/crabOff3015_C.root"); + + + TCanvas *c1 = new TCanvas("c1","c1",800,500); + c1->cd(1); + + // 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)"); + + // define cut(s) +// const char* eventcut="event%3==0 &&"; + const char* eventcut=""; + const char varcut[512]; + // sprintf(varcut,"size>%f && size<%f && length>%f &&length<%f && width>%f && width<%f",cutsize,upcutsize,lengthlowcut,lengthcut,widthlowcut,widthcut); + sprintf(varcut,"size>%f && size<%f && dist>0.2 && dist<1.1 && length>%f &&length<%f && width>%f && width<%f",cutsize,upcutsize,lengthlowcut,lengthcut,widthlowcut,widthcut); +// sprintf(varcut,"size>%f && dist>0.2 && dist<0.8 && length<0.21 && width<0.125",cutsize,upcutsize,lengthcut,widthcut); +// sprintf(varcut,"size>%f && size<%f && dist>0.2 && dist<0.8 && length<0.21 && width<0.125",cutsize,upcutsize); +// sprintf(varcut,"size>%f && size<%f && dist>0.2 && dist<0.8 && length<0.21",cutsize,upcutsize); +// sprintf(varcut,"size>%f && size<%f && dist>0.2 && dist<0.8",cutsize,upcutsize); +// sprintf(varcut,"size>%f && size<%f",cutsize,upcutsize) +// sprintf(varcut,""); + Char_t cut[256]; + Char_t cutoff[256]; + strcat(cut,eventcut); + strcat(cut,varcut); + strcat(cutoff,varcut); + + // fill on/off histos + TH1F* onhisto = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90); + TH1F* offhisto = new TH1F("OffHist","Alpha Plot",nbins,0,90); + ton->Draw("alpha>>OnHist",cut); + toff->Draw("alpha>>OffHist",cutoff); + cout << "Applied cut: " << cut << endl; + + // 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 << "Excess: " << sig-bg << endl; + cout << "N bkg: " << bg << endl; + 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(); + + // produce ps + // char psname[50]; + // sprintf(psname,"%s.ps",var); + c1->Print(psname); +} + + Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/plotOptimal.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/plotOptimal.C (revision 4056) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/plotOptimal.C (revision 4056) @@ -0,0 +1,115 @@ +//////////////////////////////////////////////////////////////////////////////////// +// +// _____ Plot Optimized cuts _____ +// +// Take as input an output from programs/optimizeCuts and plot the results +// You only need to specify the input file name +// +// Jose Flix +// Javier Rico +//////////////////////////////////////////////////////////////////////////////////// + + +void plotOptimal(TString fname="optimize.out") +{ + gROOT->Reset(); + gStyle->SetCanvasColor(0); + gStyle->SetCanvasBorderMode(0); + gStyle->SetPadBorderMode(0); + gStyle->SetFrameBorderMode(0); + gStyle->SetOptStat(000000); + gStyle->SetPalette(1); + + ifstream ifun(fname.Data()); + Float_t width,length,excess,significance; + + Float_t wmin=9999999; + Float_t lmin=9999999; + Float_t wmax=-9999999; + Float_t lmax=-9999999; + Float_t emax=-9999999; + Float_t smax=-9999999; + + Float_t dw=-1; + Float_t dl=-1; + Float_t prevw=-1; + Float_t prevl=-1; + + // first loop to get parameters + Int_t i=0; + while(ifun >> width) + { + // read data + ifun >> length; + ifun >> excess; + ifun >> significance; + + // keep max/min + if(widthwmax) wmax=width; + if(length>lmax) lmax=length; + if(significance>smax) + { + smax=significance; + emax=excess; + } + + if(width!=prevw) dw=width-prevw; + if(length!=prevl) dl=length-prevl; + + prevw=width; + prevl=length; + i++; + } + ifun.close(); + + cout << "wmin: " << wmin << ", wmax: " << wmax+dw << ", lmin: " << lmin << ", lmax: " << lmax+dl << endl; + + // define histos + Int_t nx = TMath::Nint((wmax-wmin)/dw)+1; + Int_t ny = TMath::Nint((lmax-lmin)/dl)+1; + TH2F *ex = new TH2F("Excess","Excess",nx,wmin,wmax+dw,ny,lmin,lmax+dl); + TH2F *sig = new TH2F("Significance","Significance",nx,wmin,wmax+dw,ny,lmin,lmax+dl); + + cout << "nx: " << nx << ", ny: " << ny << endl; + cout << "dw: " << dw << ", dl: " << dl << endl; + + ifstream ifun2; + ifun2.open(fname.Data()); + + // second loop to fill the histograms + while(ifun2 >> width) + { + // read data + ifun2 >> length; + ifun2 >> excess; + ifun2 >> significance; + + cout << "Width: " << width << '\t' << "Length: " << length << '\t' << "Excess: " << excess << '\t' << "Significance: " << significance; + if(significance==smax) cout << " *****"; + cout << endl; + + Int_t binx = TMath::Nint(((width-wmin)/dw)); + Int_t biny = TMath::Nint(((length-lmin)/dl)); + + // cout << "binnx: " << binx+1 << ", biny: " << biny+1 << endl; + + ex->SetBinContent(binx+1,biny+1,excess+.00001); + sig->SetBinContent(binx+1,biny+1,significance+.00001); + } + + // plot + TCanvas *c3 = new TCanvas("c3","c3",0,0,500,750); + c3->Divide(1,2); + + c3->cd(1); + ex->GetXaxis()->SetTitle("Width (degrees)"); + ex->GetYaxis()->SetTitle("Length (degrees)"); + ex->Draw("colz"); + + c3->cd(2); + sig->GetXaxis()->SetTitle("Width (degrees)"); + sig->GetYaxis()->SetTitle("Length (degrees)"); + sig->Draw("colz"); +} Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/makeHillas.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/programs/makeHillas.cc (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mifae/programs/makeHillas.cc (revision 4056) @@ -141,6 +141,4 @@ return; - // MPedestalCam pedcam = pedloop.GetPedestalCam(); - /*****************************/ /* SECOND LOOP: CALIBRATION */ @@ -155,6 +153,4 @@ if (!calloop.Process(pedloop.GetPedestalCam())) return; - - // MCalibrationChargeCam calcam = calloop.GetCalibrationCam(); /************************************************************************/ Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizeCuts.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizeCuts.cc (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizeCuts.cc (revision 4056) @@ -131,5 +131,5 @@ 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; + cout << "General cut " << gencut << " (" << evecut << " for on-data)" << endl; // loop on widht and length cuts Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizecuts.datacard =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizecuts.datacard (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mifae/programs/optimizecuts.datacard (revision 4056) @@ -7,21 +7,21 @@ // Input file name pattern (On data) -ONFILES /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OnRotateB.root +ONFILES /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OnRotateNoCalB.root // Input file name pattern (Off data) -OFFFILES /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OffRotateA-C.root +OFFFILES /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OffRotateNoCalA-C.root // output file name -OUTFILE ./optimizePrueba.root +OUTFILE ./optimize_fine.out -// Previous cuts (size in size units, distance in deg) -SIZECUT 1800 +// Preliminar cuts (size in size units, distance in deg) +SIZECUT 2000 DISTCUTS 0.2 1.1 // Length initial, final and step values -LENGTHCUTS 0.01 0.3 0.005 +LENGTHCUTS 0.1 0.4 0.005 // Width initial, final and step values -WIDTHCUTS 0.01 0.3 0.005 +WIDTHCUTS 0.0 0.3 0.005 Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/srcPos.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/programs/srcPos.cc (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mifae/programs/srcPos.cc (revision 4056) @@ -44,5 +44,5 @@ Float_t xsrcpos=0.; Float_t ysrcpos=0.; -Float_t mjdpos=53050; +Double_t mjdpos=53050; Bool_t kRotate=1; Bool_t kSrcPolicy=kTRUE; Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/srcpos.datacard =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/programs/srcpos.datacard (revision 4055) +++ trunk/MagicSoft/Mars/mtemp/mifae/programs/srcpos.datacard (revision 4056) @@ -1,8 +1,8 @@ // Maximun number of on and off events to be processed) -NEVENTS 9999999 +NEVENTS 20000 // Input file name pattern -INPUTFILES /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OnB.root +INPUTFILES /mnt/users/jrico/magic/mars/Mars_Standard02/mtemp/mifae/hillas/mrk20040215OnNoCalB.root // output file name