Index: /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc (revision 5312) +++ /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc (revision 5313) @@ -115,4 +115,5 @@ #include "MSkyPlot.h" +#include #include #include @@ -190,4 +191,9 @@ fHistOn.SetYTitle("y [\\circ]"); + fHistSignifGaus.SetName("SignifDistrib"); + fHistSignifGaus.SetTitle("Distribution of the significances from the sky plot"); + fHistSignifGaus.SetXTitle("significance"); + fHistSignifGaus.SetYTitle("counts"); + // set some values for the sky plot geometry: fMinXGrid =-1.; // [deg] @@ -274,8 +280,13 @@ fDistLo[7] = 0.0; } + + fAlphaHName = "alphaplot.root"; + fSkyHName = "skyplot.root"; } void MSkyPlot::ReadCuts(const TString parSCinit="OptSCParametersONOFFThetaRange0_1570mRad.root") { + +cout << " parameters read from file: " << parSCinit << endl; //-------------------------------- // create container for the supercut parameters @@ -390,5 +401,5 @@ fRepDrive = (MReportDrive*)plist->FindObject(AddSerialNumber("MReportDrive")); if (!fRepDrive) - *fLog << warn << "MReportDrive not found... no sky map." << endl; + *fLog << warn << "MReportDrive not found... could be problem for sky map." << endl; @@ -397,7 +408,9 @@ *fLog << warn << "MSrcPosCam not found... no sky map." << endl; +/* fPntPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MPntPosCam")); if (!fPntPosCam) *fLog << warn << "MPntPosCam not found... no sky map." << endl; +*/ fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos")); @@ -407,5 +420,5 @@ fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime")); if (!fTime) - *fLog << warn << "MTime not found... no sky map." << endl; + *fLog << warn << "MTime not found... could be problem for sky map." << endl; fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory")); @@ -443,32 +456,35 @@ fNumalphahist = (int) (fNumStepsX * fNumStepsY + 0.5); - fHistSignif.SetName("SPSignif2ndOrder"); - fHistSignif.SetTitle("Sky Plot of significance (2nd order fit)"); + // fHistSignif.SetName("SPSignif2ndOrder"); + // fHistSignif.SetTitle("Sky Plot of significance (2nd order fit)"); fHistSignif.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid); - fHistSignif.SetDirectory(NULL); + // fHistSignif.SetDirectory(NULL); fHistSignif.SetFillStyle(4000); fHistSignif.UseCurrentStyle(); - fHistNexcess.SetName("SPNex2ndOrder"); - fHistNexcess.SetTitle("Sky Plot of Number of excess events (2nd order fit)"); + // fHistNexcess.SetName("SPNex2ndOrder"); + // fHistNexcess.SetTitle("Sky Plot of Number of excess events (2nd order fit)"); fHistNexcess.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid); - fHistNexcess.SetDirectory(NULL); + // fHistNexcess.SetDirectory(NULL); fHistNexcess.SetFillStyle(4000); fHistNexcess.UseCurrentStyle(); - fHistOn.SetName("SPOnCounted"); - fHistOn.SetTitle("Sky Plot of ON events (counted)"); + // fHistOn.SetName("SPOnCounted"); + // fHistOn.SetTitle("Sky Plot of ON events (counted)"); fHistOn.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid); - fHistOn.SetDirectory(NULL); + // fHistOn.SetDirectory(NULL); fHistOn.SetFillStyle(4000); fHistOn.UseCurrentStyle(); + //fHistSignifGaus.SetName("SignifDistrib"); + fHistSignifGaus.SetTitle("Distribution of the significances from the sky plot"); + fHistSignifGaus.SetBins(100, -10., 10.); // prepare alpha plots // temp histogram for an alpha plot - TH1I *temp = new TH1I("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge); + TH1D *temp = new TH1D("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge); temp->SetDirectory(NULL); fHistAlphaBinWidth = temp->GetBinWidth(1); @@ -499,6 +515,16 @@ if (fSetCenter==kTRUE) { - fRa0 = fRepDrive->GetRa(); // [h] - fDec0 = fRepDrive->GetDec(); // [deg] + if (fRepDrive) + { + fRa0 = fRepDrive->GetRa(); // [h] + fDec0 = fRepDrive->GetDec(); // [deg] + if (fRa0 < 0. || fRa0 > 24. || fDec0 < -90. || fDec0 > 90. || (fRa0==0 && fDec0==0)) return kTRUE; // temp!, should be changed + } + else + { + fRa0 = fPointPos->GetRa(); + fDec0 = fPointPos->GetDec(); + } + *fLog << inf << "Ra (center) = " << fRa0 << ", Dec = " << fDec0 << endl; fSetCenter=kFALSE; } @@ -518,11 +544,30 @@ Double_t rho = 0; if (fTime && fObservatory && fPointPos) -{ + { rho = fPointPos->RotationAngle(*fObservatory, *fTime); -*fLog << inf << " rho = " << rho << ", Zd = " << fPointPos->GetZd() << - ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl; -} - //if (fPointPos) - // rho = fPointPos->RotationAngle(*fObservatory); +//*fLog << inf << " rho = " << rho*180./TMath::Pi() << ", Zd = " << fPointPos->GetZd() << +// ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl; + + // => coord. system: xtilde, ytilde with center in RA_0, DEC_0 + +// Get Rot. Angle: + cosgam = TMath::Cos(rho); + singam = TMath::Sin(rho); +// Get x_0, y_0 for RA_0,DEC_0 of the current event + } + else + { +// rho = fPointPos->RotationAngle(*fObservatory); + Double_t theta, phi, sing, cosg; + theta = fPointPos->GetZd()*TMath::Pi()/180.; + phi = fPointPos->GetAz()*TMath::Pi()/180.; +// printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415); + fObservatory->RotationAngle(theta, phi, sing, cosg); + cosgam = cosg; + singam = sing; + } + // if (fPointPos) + // rho = fPointPos->RotationAngle(*fObservatory); + /* //TEMP @@ -556,12 +601,10 @@ */ - // => coord. system: xtilde, ytilde with center in RA_0, DEC_0 - -// Get Rot. Angle: - cosgam = TMath::Cos(rho); - singam = TMath::Sin(rho); -// Get x_0, y_0 for RA_0,DEC_0 of the current event - x_0 = fPntPosCam->GetX()*fMm2Deg; - y_0 = fPntPosCam->GetY()*fMm2Deg; +/* + x_0 = fPntPosCam->GetX()*fMm2Deg; + y_0 = fPntPosCam->GetY()*fMm2Deg; +*/ + x_0 = 0.; + y_0 = 0.; Int_t index = 0; // index for alpha histograms @@ -598,5 +641,5 @@ lgsize = logsize-log3000; lgsize2 = lgsize*lgsize; - dist2 = (distpar-fDistOffset)*(distpar-fDistOffset); + dist2 = distpar*distpar - fDistOffset*fDistOffset; if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) && @@ -609,7 +652,9 @@ { // gamma candidates! +//*fLog << "Length : " << fHillas->GetLength()*fMm2Deg << ", Width : " << fHillas->GetWidth()*fMm2Deg << endl; +//*fLog << "distpar: " << distpar << ", Size : " << fHillas->GetSize() << endl; fHistAlpha[index].Fill(alphapar); - index++; } + index++; } } @@ -628,4 +673,7 @@ const Int_t onbinsalpha = TMath::Nint(fAlphaONMax/fHistAlphaBinWidth); +// fit with gaus + fHistSignifGaus.Fit("gaus"); + // fit function for the background TF1 * fitbgpar = new TF1("fbgpar", "[0]*x*x + [1]", fAlphaBgLow, fAlphaBgUp); @@ -635,5 +683,5 @@ numdegfreed = (fAlphaBgUp - fAlphaBgLow) / fHistAlphaBinWidth - 2.; - vector ::iterator alpha_iterator; + vector ::iterator alpha_iterator; int index2 = 0; // index of the TH2F histograms @@ -643,5 +691,6 @@ // find the bin in the 2dim histogram nrow = index2/fHistOn.GetNbinsX() + 1; // row of the histogram (y) - ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1; // column of the histogram (x) + ncolumn = index2%fHistOn.GetNbinsX()+1; // column of the histogram (x) + //ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1; // column of the histogram (x) // number of ON events below fAlphaONMax @@ -659,5 +708,5 @@ // estimate Noff from the fit: - Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * pow(fAlphaONMax,3.) + + Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * TMath::Power(fAlphaONMax,3.) + (fitbgpar->GetParameter(1)) * fAlphaONMax) / fHistAlphaBinWidth; @@ -671,4 +720,5 @@ fHistSignif.SetBinContent(ncolumn, nrow, Sign); // fill in the fHistSignif + fHistSignifGaus.Fill(Sign); alpha_iterator++; @@ -677,4 +727,5 @@ //temp +/* Double_t decl, hang; MStarCamTrans mstarc(*fGeomCam, *fObservatory); @@ -683,11 +734,15 @@ *fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl; *fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl; +*/ //endtemp + // save alpha plots: - if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots(); +// TString stri1 = "alphaplots.root"; + if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots(fAlphaHName); // save sky plots: - if(kSaveSkyPlots==kTRUE) SaveSkyPlots(); +// TString stri2 = "skyplots.root"; + if(kSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName); fHistAlpha.clear(); @@ -721,5 +776,5 @@ -void MSkyPlot::SaveSkyPlots(const TString skyplotfilename="skyplots.root") +void MSkyPlot::SaveSkyPlots(TString skyplotfilename) { TFile rootfile(skyplotfilename, "RECREATE", @@ -728,4 +783,5 @@ fHistNexcess.Write(); fHistOn.Write(); + fHistSignif.Write(); // from Wolfgang: @@ -733,17 +789,22 @@ // Dec-Hour lines Double_t rho, sinrho, cosrho; -// if (fTime && fObservatory && fPointPos) + Double_t theta, phi, sing, cosg; +// do I need it? + if (fTime && fObservatory && fPointPos) + { rho = fPointPos->RotationAngle(*fObservatory, *fTime); - sinrho=TMath::Sin(rho); - cosrho=TMath::Cos(rho); - - Double_t theta, phi, sing, cosg; - theta = fPointPos->GetZd()*TMath::Pi()/180.; - phi = fPointPos->GetAz()*TMath::Pi()/180.; + sinrho=TMath::Sin(rho); + cosrho=TMath::Cos(rho); + } + + theta = fPointPos->GetZd()*TMath::Pi()/180.; + phi = fPointPos->GetAz()*TMath::Pi()/180.; // printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415); - fObservatory->RotationAngle(theta, phi, sing, cosg); + fObservatory->RotationAngle(theta, phi, sing, cosg); *fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl; *fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl; + sinrho=sing; + cosrho=cosg; Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0; // [mm] @@ -756,7 +817,8 @@ Double_t declin, hangle; // [deg], [h] MStarCamTrans mstarc(*fGeomCam, *fObservatory); - mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle); - - TGraph *graph1; + if (fRepDrive) mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle); + else mstarc.LocToCel(theta*180./TMath::Pi(),phi*180./TMath::Pi(),declin, hangle); // NOT GOOD! + +// std::vector graphvec; TLatex *pix; @@ -776,10 +838,22 @@ // direction for the center of the camera Double_t dec0 = declin; - Double_t h0 = hangle*360./24.; + Double_t h0 = hangle*360./24.; //deg + Double_t RaHOffset = fRepDrive->GetRa() - h0; //trans.LocToCel(theta0, phi0, dec0, h0); + gStyle->SetOptFit(0); + gStyle->SetOptStat(0); + gStyle->SetPalette(1); TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600); + c1->Divide(2,2); c1->cd(1); + fHistSignif.Draw("colz"); + c1->cd(2); + fHistNexcess.Draw("colz"); + c1->cd(3); fHistOn.Draw("colz"); + c1->cd(4); + gPad->SetLogy(); + fHistSignifGaus.Draw(); //----- lines for fixed dec ------------------------------------ @@ -796,5 +870,5 @@ dh = 360.0/((Double_t)(Nh-1)); - +// start to copy for (Int_t j=0; jcd(2); - graph1 = new TGraph(Nh,xh,yh); - graph1->SetLineColor(j+1); + TGraph * graph1 = new TGraph(Nh,xh,yh); + //graph1->SetLineColor(j+1); + graph1->SetLineColor(36); graph1->SetLineStyle(1); - graph1->SetMarkerColor(j+1); + graph1->SetLineWidth(1); + //graph1->SetMarkerColor(j+1); + graph1->SetMarkerColor(1); graph1->SetMarkerSize(.2); graph1->SetMarkerStyle(20); - graph1->Draw("PL"); + + c1->cd(1); + graph1->Draw("L"); + c1->cd(2); + graph1->Draw("L"); + c1->cd(3); + graph1->Draw("L"); //delete graph1; +// graphvec.push_back(*graph1); +// graphvec[j].Draw("L"); sprintf(tit,"Dec = %6.2f", dec[j]); @@ -852,11 +935,11 @@ ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)j) *(yup-ylo)/20.0; pix = new TLatex(xtxt, ytxt, tit); - pix->SetTextColor(j+1); + pix->SetTextColor(36); pix->SetTextSize(.03); - pix->Draw(""); + //pix->Draw(""); //delete pix; } - +//stop copy //----- lines for fixed hour angle ------------------------------------ @@ -872,5 +955,5 @@ dh1 = 360.0/((Double_t)Nh1); - +// start copy for (Int_t j=0; jcd(2); - graph1 = new TGraph(count,xd,yd); - graph1->SetLineColor(k+1); - graph1->SetLineStyle(2); - graph1->SetMarkerColor(k+1); + TGraph * graph1 = new TGraph(count,xd,yd); + //graph1->SetLineColor(k+1); + graph1->SetLineColor(36); + graph1->SetLineWidth(2); + graph1->SetLineStyle(1); + //graph1->SetMarkerColor(k+1); + graph1->SetMarkerColor(1); graph1->SetMarkerSize(.2); graph1->SetMarkerStyle(20); - graph1->Draw("PL"); - //delete graph1; - - sprintf(tit,"H = %6.2f", h1[k]); + c1->cd(1); + graph1->Draw("L"); + c1->cd(2); + graph1->Draw("L"); + c1->cd(3); + graph1->Draw("L"); + + sprintf(tit,"RA = %6.2f", fRa0 + (h0 - h1[k])); Double_t xtxt = xlo + (xup-xlo)*0.75; Double_t ytxt = ylo + (yup-ylo)*0.80 - ((Double_t)k) *(yup-ylo)/20.0; pix = new TLatex(xtxt, ytxt, tit); - pix->SetTextColor(k+1); + pix->SetTextColor(36); pix->SetTextSize(.03); - pix->Draw(""); + //pix->Draw(""); //delete pix; @@ -936,22 +1026,23 @@ // c1->cd(2); - sprintf(tit,"Dec0 = %6.2f H0 = %6.2f", dec0, h0); - xtxt = xlo + (xup-xlo)*0.05; - ytxt = ylo + (yup-ylo)*0.92; + sprintf(tit,"Dec0 = %6.2f [deg] Ra0 = %6.2f [h]", dec0, fRa0); + xtxt = xlo + (xup-xlo)*0.05 + 0.80; + ytxt = ylo + (yup-ylo)*0.75; pix = new TLatex(xtxt, ytxt, tit); pix->SetTextColor(1); - pix->SetTextSize(.06); + pix->SetTextSize(.05); + c1->cd(1); + pix->Draw(""); + c1->cd(2); + pix->Draw(""); + c1->cd(3); pix->Draw(""); //delete pix; - sprintf(tit," [deg] [deg]"); - xtxt = xlo + (xup-xlo)*0.05; - ytxt = ylo + (yup-ylo)*0.86; - pix = new TLatex(xtxt, ytxt, tit); - pix->SetTextColor(1); - pix->SetTextSize(.06); - pix->Draw(""); - c1->Write(); + // we suppose that the {skyplotfilename} ends with .root + Int_t sizeofout = skyplotfilename.Sizeof(); + TString outps = skyplotfilename.Remove(sizeofout-5,5) + "ps"; + c1->Print(outps); // temporary!!!!! TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600); @@ -967,14 +1058,15 @@ delete c1; delete c2; - delete graph1; delete pix; + + } -void MSkyPlot::SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root") +void MSkyPlot::SaveAlphaPlots(const TString alphaplotfilename) { TFile rootfile(alphaplotfilename, "RECREATE", "all the alpha plots"); - vector ::iterator alpha_iterator; + vector ::iterator alpha_iterator; int index = 0; // index of the TH2F histograms Char_t strtitle[100]; @@ -989,5 +1081,6 @@ nrow = index/fHistOn.GetNbinsX() + 1; // row of the histogram (y) - ncolumn = TMath::Nint(fmod(index,fHistOn.GetNbinsX()))+1; // column of the histogram (x) + //ncolumn = TMath::Nint(fmod(index,fHistOn.GetNbinsX()))+1; // column of the histogram (x) + ncolumn = index%fHistOn.GetNbinsX()+1; // column of the histogram (x) xpos = fMinXGrid + fBinStepGrid*(ncolumn-1); ypos = fMinYGrid + fBinStepGrid*(nrow-1); Index: /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h (revision 5312) +++ /trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h (revision 5313) @@ -17,5 +17,5 @@ class TH2D; -class TH1I; +class TH1D; class MParList; @@ -59,5 +59,5 @@ // Float_t fMaxLD; // Maximum distance in percent of dist - vector fHistAlpha; // vector of histograms for alpha + std::vector fHistAlpha; // vector of histograms for alpha Int_t fNumalphahist; // number of histograms for alpha Int_t fNumBinsAlpha; @@ -74,4 +74,5 @@ TH2D fHistNexcess; // sky plot of number of excess events vs. x and y TH2D fHistOn; // sky plot of events below fAlphaONMax vs. x and y + TH1D fHistSignifGaus; // distribution of significance Bool_t fSetCenter; // used to set the center of these histograms once Double_t fRa0; @@ -85,6 +86,6 @@ Float_t fMaxYGrid; // [deg] , lower edge of the skyplot Float_t fBinStepGrid; // [deg] , grid size - Float_t fNumStepsX; // number of bins in x - Float_t fNumStepsY; // number of bins in y + Int_t fNumStepsX; // number of bins in x + Int_t fNumStepsY; // number of bins in y @@ -101,4 +102,7 @@ Double_t fRa; Double_t fDec; + + TString fAlphaHName; // name for histogram with alpha plots + TString fSkyHName; // name for histogram with sky plots Int_t DistancetoPrimitive(Int_t px, Int_t py); @@ -118,4 +122,6 @@ void SetMinDist(Float_t dist) { fMinDist = dist; } // Absolute minimum distance void SetMaxDist(Float_t dist) { fMaxDist = dist; } // Absolute maximum distance + void SetSizeMin(Float_t size) { fSizeMin = size; } // Absolute minimum Size + void SetSizeMax(Float_t size) { fSizeMax = size; } // Absolute maximum Size // void SetMinLD(Float_t ratio) { fMinLD = ratio; } // Minimum ratio between length/dist // void SetMaxLD(Float_t ratio) { fMaxLD = ratio; } // Maximum ratio between length/dist @@ -123,9 +129,19 @@ void SetSkyPlot(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t step); Double_t CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2); - void SaveSkyPlots(const TString skyplotfilename="skyplots.root"); - void SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root"); + //void SaveSkyPlots(const TString skyplotfilename="skyplots.root"); + + void SetOutputAlphaName(TString outname1) { fAlphaHName = outname1; } + void SaveSkyPlots(TString stri); + + + //void SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root"); + void SetOutputSkyName(TString outname2) { fSkyHName = outname2; } + void SaveAlphaPlots(const TString stri2); void SetAlphaCut(Float_t alpha); + //std::vector ::iterator GetAlphaPlots() { std::vector ::iterator iter; return iter = fHistAlpha.begin() ; } + std::vector GetAlphaPlots() { return fHistAlpha ; } + void Draw(Option_t *option=""); Index: /trunk/MagicSoft/Mars/mtemp/mmpi/macros/MySkyPlot.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mmpi/macros/MySkyPlot.C (revision 5313) +++ /trunk/MagicSoft/Mars/mtemp/mmpi/macros/MySkyPlot.C (revision 5313) @@ -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 26.09.2004 mazin@mppmu.mpg.de +//********************************************************** +// this macro : +// 1) take Hillas Parameters and produce a sky plot +//********************************************************** + +gROOT->Reset(); + +//************************************************************************ +void MySkyPlot(const char *onfile) +{ + + TString fileON = onfile; + + // name of input root file + TString filenamein(fileON); + + + // set up the task list: + //----------------------------------------------------------- + MTaskList tlist; + MParList plist; + + + TString XPosMaxString ("MSrcPosCam.fX<0.1"); + MContinue XPosMaxCut(XPosMaxString); + XPosMaxCut.SetInverted(); + + TString XPosMinString ("MSrcPosCam.fX>-0.1"); + MContinue XPosMinCut(XPosMinString); + XPosMinCut.SetInverted(); + + TString YPosMaxString ("MSrcPosCam.fY<0.1"); + MContinue YPosMaxCut(YPosMaxString); + YPosMaxCut.SetInverted(); + + TString YPosMinString ("MSrcPosCam.fY>-0.1"); + MContinue YPosMinCut(YPosMinString); + YPosMinCut.SetInverted(); + + + + MSkyPlot skyplot; + skyplot.SetMinDist(0.2); + //skyplot.ReadCuts("/local/supercuts/OptSCParametersONOFFThetaRange0_1570mRad.root"); + //skyplot.ReadCuts("/remote/home/pc1/operator/mazin/supercuts/2004_09_13/Zd20_60all/Cuts.root"); + //skyplot.ReadCuts("/local/supercuts/2004_09_13/OptDynCuts2000CrabNebula.root"); + //skyplot.ReadCuts("/remote/home/pc17/operator/OnlineA/supercuts/OptDynCuts2000CrabNebula.root"); + //skyplot.ReadCuts("/home/pcmagic04/pratik/crab_analysis/sept2k4/gt2000/RootFileDynCuts_crabsept_GREEN_15o30deg.root"); + skyplot.ReadCuts("/home/pcmagic04/pratik/crab_analysis/sept2k4/gt2000/RootFileDynCuts_20to30deg_alldata.root"); + skyplot.SetSkyPlot(-1.,1.,-1.,1.,0.05); + skyplot.SetSizeMin(2000.); + skyplot.SetSizeMax(1e6); + skyplot.SetOutputAlphaName("Crab21_UV.root"); + skyplot.SetOutputSkyName("Crab21_UV.root"); + + + MObservatory mobs; + + plist.AddToList(&mobs); + + // geometry is needed in MHHillas... classes + const MGeomCam *fGeom = + (MGeomCam*)plist->FindCreateObj("MGeomCamMagic", "MGeomCam"); + + //------------------------------------------- + // create the tasks which should be executed + // + + MReadReports read; + read.AddTree("Drive"); +// read.AddTree("CC"); +// read.AddTree("Events"); + read.AddTree("Events","MTime.",kTRUE); + + read.AddFile(filenamein); + + read.AddToBranchList("MReportDrive.*"); +// read.AddToBranchList("MRawEvtHeader.*"); + + + //MReadMarsFile read("Events", filenamein); + //MReadMarsFile read(filenamein); + //read.DisableAutoScheme(); + + MGeomApply apply; + + + + //***************************** + // entries in MParList + + plist.AddToList(&tlist); + +// pliston.AddToList(&source); + + //***************************** + // entries in MTaskList + + + tlist.AddToList(&read); + tlist.AddToList(&apply); + + tlist.AddToList(&XPosMaxCut); + tlist.AddToList(&XPosMinCut); + + tlist.AddToList(&YPosMaxCut); + tlist.AddToList(&YPosMinCut); + + tlist.AddToList(&skyplot); + + //***************************** + + //------------------------------------------- + // Execute event loop + // + MProgressBar bar; + MEvtLoop evtloop; + evtloop.SetParList(&plist); + evtloop.SetProgressBar(&bar); +/* + Int_t maxevents = -1; + //Int_t maxevents = 1000; + if ( !evtloop.Eventloop(maxevents) ) + return; +*/ + if (!evtloop.PreProcess()) + return; + + Int_t counter=0; + while(tlist.Process()) + { + counter++; + if(counter%10000 == 0) cout << "events = " << counter << endl; +// MReportDrive * mrepdrive = (MReportDrive *)pliston.FindObject("MReportDrive"); +// cout << "Ra " << mrepdrive->GetRa() << ", Dec " << mrepdrive->GetDec() << endl; +// MHillasSrc * mhs = (MHillasSrc*)pliston.FindObject("MHillasSrc"); +// mhs->Print(); + } + + evtloop.PostProcess(); + + tlist.PrintStatistics(0, kTRUE); + +} + + +