Index: trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc (revision 6189) +++ trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc (revision 6190) @@ -115,5 +115,4 @@ #include "MSkyPlot.h" -#include #include #include @@ -135,4 +134,5 @@ #include "MHillasExt.h" #include "MNewImagePar.h" +#include "MHadronness.h" #include "MTime.h" #include "MObservatory.h" @@ -151,4 +151,6 @@ #include "MLogManip.h" +#include + ClassImp(MSkyPlot); @@ -160,9 +162,19 @@ // MSkyPlot::MSkyPlot(const char *name, const char *title) - : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1) -// fAlphaCut(12.5), BgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1) -{ - -*fLog << warn << "entering default constructor in MSkyPlot" << endl; + : fGeomCam(NULL), + fTime(NULL), + fPointPos(NULL), + fRepDrive(NULL), + fSrcPosCam(NULL), + fPntPosCam(NULL), + fObservatory(NULL), + fHillas(NULL), + fHillasExt(NULL), + fHillasSrc(NULL), + fNewImagePar(NULL), + fMm2Deg(-1) +{ + + *fLog << warn << "entering default constructor in MSkyPlot" << endl; // // set the name and title of this object @@ -176,4 +188,10 @@ fHistNexcess.SetDirectory(NULL); fHistOn.SetDirectory(NULL); + fHistSignifGaus.SetDirectory(NULL); + + fHistSignif.UseCurrentStyle(); + fHistNexcess.UseCurrentStyle(); + fHistOn.UseCurrentStyle(); + fHistSignifGaus.UseCurrentStyle(); fHistSignif.SetName("SkyPlotSignif"); @@ -210,4 +228,6 @@ fGridFineBin = 0.01; // degrees +// elapsed time: + fElaTime = 0.; // some filter cuts: @@ -217,4 +237,5 @@ fMaxDist = 1.4; // dist max cut (ever possible) fMinDist = 0.1; // dist max cut (ever possible) + fHadrCut = 0.2; // hadronness cut fNumBinsAlpha = 36; // number of bins for alpha histograms @@ -282,11 +303,48 @@ } - kSaveAlphaPlots=kTRUE; - kSaveSkyPlots=kTRUE; - kSaveNexPlot=kTRUE; + fSaveAlphaPlots=kTRUE; + fSaveSkyPlots=kTRUE; + fSaveNexPlot=kTRUE; + fUseRF=kFALSE; fAlphaHName = "alphaplot.root"; fSkyHName = "skyplot.root"; fNexHName = "Nexcess.gif"; -} + + fHistAlpha = new TOrdCollection(); + fHistAlpha->SetOwner(); + +} + +MSkyPlot::~MSkyPlot() +{ + + if (fHistAlpha) + delete fHistAlpha; +} + +// -------------------------------------------------------------------------- +// +// Get i-th hist +// +TH1D *MSkyPlot::GetAlphaPlot(Int_t i) +{ + if (GetSize() == 0) + return NULL; + + return static_cast(i==-1 ? fHistAlpha->At(GetSize()/2+1) : fHistAlpha->At(i)); +} + +// -------------------------------------------------------------------------- +// +// Get i-th hist +// +const TH1D *MSkyPlot::GetAlphaPlot(Int_t i) const +{ + if (GetSize() == 0) + return NULL; + + return static_cast(i==-1 ? fHistAlpha->At(GetSize()/2+1) : fHistAlpha->At(i)); +} + void MSkyPlot::ReadCuts(const TString parSCinit="OptSCParametersONOFFThetaRange0_1570mRad.root") @@ -412,6 +470,4 @@ } - - // -------------------------------------------------------------------------- // @@ -423,5 +479,7 @@ Int_t MSkyPlot::PreProcess(MParList *plist) { -*fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl; + + *fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl; + fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam"); if (!fGeomCam) @@ -438,5 +496,5 @@ - fSrcPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam")); + fSrcPosCam = (MSrcPosCam*)plist->FindCreateObj(AddSerialNumber("MSrcPosCam")); if (!fSrcPosCam) *fLog << warn << "MSrcPosCam not found... no sky map." << endl; @@ -477,5 +535,12 @@ *fLog << err << "MNewImagePar not found... no sky map." << endl; - + if(fUseRF) + { + fHadron = (MHadronness*)plist->FindObject(AddSerialNumber("MHadronness")); + if (!fHadron) + *fLog << err << "MHadronness not found although specified... no sky map." << endl; + + *fLog << inf << "Hadronness cut set to : " << fHadrCut << endl; + } // FIXME: Because the pointing position could change we must check @@ -495,7 +560,5 @@ fHistSignif.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid); - // fHistSignif.SetDirectory(NULL); fHistSignif.SetFillStyle(4000); - fHistSignif.UseCurrentStyle(); // fHistNexcess.SetName("SPNex2ndOrder"); @@ -503,7 +566,5 @@ fHistNexcess.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid); - // fHistNexcess.SetDirectory(NULL); fHistNexcess.SetFillStyle(4000); - fHistNexcess.UseCurrentStyle(); // fHistOn.SetName("SPOnCounted"); @@ -511,7 +572,5 @@ fHistOn.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid); - // fHistOn.SetDirectory(NULL); fHistOn.SetFillStyle(4000); - fHistOn.UseCurrentStyle(); //fHistSignifGaus.SetName("SignifDistrib"); @@ -519,17 +578,16 @@ fHistSignifGaus.SetBins(100, -10., 10.); - // prepare alpha plots + // create alpha histograms + for (Int_t i=0; i< fNumalphahist; i++) + { // temp histogram for an alpha plot - TH1D *temp = new TH1D("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge); - temp->SetDirectory(NULL); - fHistAlphaBinWidth = temp->GetBinWidth(1); - // create alpha histograms - for (Int_t i=0; i< fNumalphahist; i++) { - fHistAlpha.push_back(*temp); - fHistAlpha[i].SetDirectory(NULL); - } -//cout << " fHistAlpha[10].GetBinContent(5) " << fHistAlpha[10].GetBinContent(5) << endl; - - delete temp; + TH1D *temp = new TH1D("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge); + temp->SetDirectory(NULL); + fHistAlpha->AddAt(temp,i); + } + + fHistAlphaBinWidth = GetAlphaPlot()->GetBinWidth(1); +//cout << " (*fHistAlpha)[10].GetBinContent(5) " << (*fHistAlpha)[10].GetBinContent(5) << endl; + return kTRUE; } @@ -541,167 +599,186 @@ Int_t MSkyPlot::Process() { -// check whether MPointingPos comtains something: + +// *fLog << err << "MPointingPos ENTERING the process" << endl; +// *fLog << err << "MPointingPos:: fUseRF " << (int)fUseRF << endl; + // check whether MPointingPos comtains something: if (TMath::Abs(fPointPos->GetRa())<1e-3 && TMath::Abs(fPointPos->GetDec())<1e-3) - { - *fLog << warn << "MPointingPos is not filled ... event skipped" << endl; - return kTRUE; - } - - // Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside) + { + *fLog << warn << "MPointingPos is not filled ... event skipped" << endl; + return kTRUE; + } + + // Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside) if (fSetCenter==kTRUE) - { - 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; - } - - -// some filter cuts: + { + 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; + } + + // some filter cuts: if ( fHillas->GetSize() > fSizeMin && fHillas->GetSize() < fSizeMax && fNewImagePar->GetLeakage1() < fLeakMax) - { - - Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam; - Double_t dx, dy, beta, tanbeta, alphapar, distpar; - Double_t logsize, lgsize, lgsize2, dist2; - const Double_t log3000 = log(3000.); - const Float_t fDistOffset = 0.9; - - // Get camera rotation angle - Double_t rho = 0; - if (fTime && fObservatory && fPointPos) - { - rho = fPointPos->RotationAngle(*fObservatory, *fTime); -//*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 -// theta = mcevt->GetTelescopeTheta(); - 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); + { -// conclusion: diffference in rho = 7 deg - *fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ", " - << phi << ", " << cosg << ", " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl; - - Double_t a1 = 0.876627; - Double_t a3 = -0.481171; - Double_t thetaTel=theta, phiTel=phi; - - Double_t denom = 1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) + - ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) * - ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) ); - Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom; - Double_t sinal = a1 * sin(phiTel) * denom; - - *fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ", " - << phiTel << ", " << cosal << ", " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl; - - - -// END TEMP -*/ - -// make the center of the plot different from the center of the camera -/* - x_0 = fPntPosCam->GetX()*fMm2Deg; - y_0 = fPntPosCam->GetY()*fMm2Deg; -*/ - x_0 = 0.; - y_0 = 0.; - - Int_t index = 0; // index for alpha histograms - // loop over xtilde - for (Int_t gridy = 0; gridy < fNumStepsY; gridy++) - { - ysource = fMinYGrid + gridy * fBinStepGrid; - // loop over ytilde - for (Int_t gridx = 0; gridx < fNumStepsX; gridx++) + Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam; + Double_t dx, dy, beta, tanbeta, alphapar, distpar; + Double_t logsize, lgsize, lgsize2, dist2, hadr; + const Double_t log3000 = log(3000.); + const Float_t fDistOffset = 0.9; + + //Get Hadronness if available: + if(fUseRF) { - - xsource = fMinXGrid + gridx * fBinStepGrid; - - /* derotation : rotate into camera coordinates */ - /* formel: (x_cam) (cos(gam) -sin(gam)) (xtilde) (x_0) - ( ) = ( ) * ( ) + ( ) - (y_cam) (sin(gam) sin(gam)) (ytilde) (y_0) -*/ - xsourcam = cosgam * xsource - singam * ysource + x_0; - ysourcam = singam * xsource + cosgam * ysource + y_0; - - - /* end derotatiom */ -// xsourcam = xsource; -// ysourcam = ysource; - - /* calculate ALPHA und DIST according to the source position */ - dx = fHillas->GetMeanX()*fMm2Deg - xsourcam; - dy = fHillas->GetMeanY()*fMm2Deg - ysourcam; - tanbeta = dy / dx ; - beta = TMath::ATan(tanbeta); - alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi(); - distpar = sqrt( dy*dy + dx*dx ); - if(alphapar > 90.) alphapar -= 180.; - if(alphapar < -90.) alphapar += 180.; - alphapar = TMath::Abs(alphapar); - -// apply cuts - logsize = log(fHillas->GetSize()); - lgsize = logsize-log3000; - lgsize2 = lgsize*lgsize; - dist2 = distpar*distpar - fDistOffset*fDistOffset; - - if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) && - (fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2)) - if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) && - (fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2)) - if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) && - distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) && - distpar < fMaxDist && distpar > fMinDist) - { -// gamma candidates! -//*fLog << "Length : " << fHillas->GetLength()*fMm2Deg << ", Width : " << fHillas->GetWidth()*fMm2Deg << endl; -//*fLog << "distpar: " << distpar << ", Size : " << fHillas->GetSize() << endl; - fHistAlpha[index].Fill(alphapar); - } - index++; + hadr=fHadron->GetHadronness(); +// cout << " will use RF " << hadr << endl; } - } - } - + // Get camera rotation angle + Double_t rho = 0; + if (fTime && fObservatory && fPointPos) + { + rho = fPointPos->RotationAngle(*fObservatory, *fTime); + //*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 + // theta = mcevt->GetTelescopeTheta(); + 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); + + // conclusion: diffference in rho = 7 deg + // *fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ", " + // << phi << ", " << cosg << ", " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl; + + Double_t a1 = 0.876627; + Double_t a3 = -0.481171; + Double_t thetaTel=theta, phiTel=phi; + + Double_t denom = 1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) + + ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) * + ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) ); + Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom; + Double_t sinal = a1 * sin(phiTel) * denom; + + // *fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ", " + // << phiTel << ", " << cosal << ", " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl; + + // END TEMP + */ + + // make the center of the plot different from the center of the camera + /* + x_0 = fPntPosCam->GetX()*fMm2Deg; + y_0 = fPntPosCam->GetY()*fMm2Deg; + */ + x_0 = 0.; + y_0 = 0.; + + Int_t index = 0; // index for alpha histograms + // loop over xtilde + for (Int_t gridy = 0; gridy < fNumStepsY; gridy++) + { + ysource = fMinYGrid + gridy * fBinStepGrid; + // loop over ytilde + for (Int_t gridx = 0; gridx < fNumStepsX; gridx++) + { + + xsource = fMinXGrid + gridx * fBinStepGrid; + + /* derotation : rotate into camera coordinates */ + /* formel: (x_cam) (cos(gam) -sin(gam)) (xtilde) (x_0) + ( ) = - ( ) * ( ) + ( ) + (y_cam) (sin(gam) cos(gam)) (ytilde) (y_0) + */ + //xsourcam = - (cosgam * xsource - singam * ysource) + x_0; + //ysourcam = - (singam * xsource + cosgam * ysource) + y_0; + + + /* end derotatiom */ + xsourcam = xsource; + ysourcam = ysource; + + /* calculate ALPHA und DIST according to the source position */ + dx = fHillas->GetMeanX()*fMm2Deg - xsourcam; + dy = fHillas->GetMeanY()*fMm2Deg - ysourcam; + tanbeta = dy / dx ; + beta = TMath::ATan(tanbeta); + alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi(); + distpar = sqrt( dy*dy + dx*dx ); + if(alphapar > 90.) alphapar -= 180.; + if(alphapar < -90.) alphapar += 180.; + alphapar = TMath::Abs(alphapar); + + if(fUseRF) + { + +// cout << " will use RF " << hadr << endl; + if(hadr fMinDist) + { + TH1D *hist = GetAlphaPlot(index); + hist->Fill(alphapar); + } + } + else + { + // apply cuts + logsize = log(fHillas->GetSize()); + lgsize = logsize-log3000; + lgsize2 = lgsize*lgsize; + dist2 = distpar*distpar - fDistOffset*fDistOffset; + + if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) && + (fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2)) + if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) && + (fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2)) + if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) && + distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) && + distpar < fMaxDist && distpar > fMinDist) + { + // gamma candidates! + //*fLog << "Length : " << fHillas->GetLength()*fMm2Deg << ", Width : " << fHillas->GetWidth()*fMm2Deg << endl; + //*fLog << "distpar: " << distpar << ", Size : " << fHillas->GetSize() << endl; + TH1D *hist = GetAlphaPlot(index); + hist->Fill(alphapar); + } + } + index++; + } + } + } return kTRUE; } @@ -724,8 +801,8 @@ numdegfreed = (fAlphaBgUp - fAlphaBgLow) / fHistAlphaBinWidth - 2.; - vector ::iterator alpha_iterator; int index2 = 0; // index of the TH2F histograms - alpha_iterator = fHistAlpha.begin(); + TOrdCollectionIter Next(fHistAlpha); + TH1D *alpha_iterator = NULL; fHistNexcess.Reset(); @@ -734,8 +811,8 @@ fHistSignifGaus.Reset(); - while( alpha_iterator != fHistAlpha.end() ) { + while ( (alpha_iterator = (TH1D*)Next())) { // find the bin in the 2dim histogram - nrow = index2/fHistOn.GetNbinsX() + 1; // row of the histogram (y) - ncolumn = index2%fHistOn.GetNbinsX()+1; // column of the histogram (x) + nrow = index2/fHistOn.GetNbinsX() + 1; // row of the histogram (y) + ncolumn = index2%fHistOn.GetNbinsX()+1; // column of the histogram (x) //ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1; // column of the histogram (x) @@ -747,58 +824,60 @@ // fit parabola to a background region - (*alpha_iterator).Fit(fitbgpar,"EQRN"); // NWR OK????????????????????????? + alpha_iterator->Fit(fitbgpar,"EQRN"); // NWR OK????????????????????????? // get Chi2 chisquarefit = fitbgpar->GetChisquare(); if (chisquarefit/numdegfreed<2. && chisquarefit/numdegfreed>0.01); else *fLog << warn << "Fit is bad, chi2/ndf = " << chisquarefit/numdegfreed << endl; - + // estimate Noff from the fit: Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * TMath::Power(fAlphaONMax,3.) + - (fitbgpar->GetParameter(1)) * fAlphaONMax) / fHistAlphaBinWidth; - + (fitbgpar->GetParameter(1)) * fAlphaONMax) / fHistAlphaBinWidth; + Nex = Non - Noff_fit; fHistNexcess.SetBinContent(ncolumn, nrow, Nex); // fill in the fHistNexcess - + if (Noff_fit<0.) Sign = 0.; -// else Sign = LiMa17(Non,Noff_fit,1.); + // else Sign = LiMa17(Non,Noff_fit,1.); else Sign = MMath::SignificanceLiMaSigned(Non, Noff_fit, 1.); - + fHistSignif.SetBinContent(ncolumn, nrow, Sign); // fill in the fHistSignif fHistSignifGaus.Fill(Sign); - - alpha_iterator++; + index2++; } - -// fit with gaus + + // fit with gaus fHistSignifGaus.Fit("gaus","N"); - - -//temp -/* - Double_t decl, hang; - MStarCamTrans mstarc(*fGeomCam, *fObservatory); - mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang); - -*fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl; -*fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl; -*/ -//endtemp - - -// save alpha plots: -// TString stri1 = "alphaplots.root"; - if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots(fAlphaHName); - -// save sky plots: -// TString stri2 = "skyplots.root"; - if(kSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName); - -// save nex plot: - if(kSaveNexPlot==kTRUE) SaveNexPlot(fNexHName); - - fHistAlpha.clear(); - return kTRUE; + + + //temp + /* + Double_t decl, hang; + MStarCamTrans mstarc(*fGeomCam, *fObservatory); + mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang); + + *fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl; + *fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl; + */ + //endtemp + + + // save alpha plots: + // TString stri1 = "alphaplots.root"; + if(fSaveAlphaPlots==kTRUE) SaveAlphaPlots(fAlphaHName); + + // save sky plots: + // TString stri2 = "skyplots.root"; + if(fSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName); + + // save nex plot: + if(fSaveNexPlot==kTRUE) SaveNexPlot(fNexHName); + + fHistAlpha->Clear(); + + delete fitbgpar; + + return kTRUE; } @@ -843,4 +922,6 @@ gStyle->SetPadLeftMargin(0.13); + Char_t timet[100]; + TH2D tmp = fHistNexcess; tmp.SetMaximum(470); @@ -861,4 +942,10 @@ gPad->SetTheta(20); fHistNexcess.Draw("lego2"); + TLatex tu(0.5,0.8,"elapsed time:"); + tu.Draw(); + sprintf(timet,"%.1f min",fElaTime); + TLatex tut(0.5,0.7,timet); + tut.Draw(); + can.Print(nexplotname); // can.Print("test.root"); @@ -867,137 +954,137 @@ void MSkyPlot::SaveSkyPlots(TString skyplotfilename) { - TFile rootfile(skyplotfilename, "RECREATE", - "sky plots in some variations"); - fHistSignif.Write(); - fHistNexcess.Write(); - fHistOn.Write(); - fHistSignif.Write(); - -// from Wolfgang: - //-------------------------------------------------------------------- - // Dec-Hour lines - Double_t rho, sinrho, cosrho; - 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); - } - - 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); - -*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] - Double_t gridbinning = fGridBinning; - Double_t gridfinebin = fGridFineBin; - Int_t numgridlines = (Int_t)(4.0/gridbinning); - Double_t aberr = 1.07; - Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ; - - Double_t declin, hangle; // [deg], [h] - MStarCamTrans mstarc(*fGeomCam, *fObservatory); - 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; - - Char_t tit[100]; - Double_t xtxt; - Double_t ytxt; - - Double_t xlo = -534.0 * mmtodeg; - Double_t xup = 534.0 * mmtodeg; - - Double_t ylo = -534.0 * mmtodeg; - Double_t yup = 534.0 * mmtodeg; - - Double_t xx, yy; - - - // direction for the center of the camera - Double_t dec0 = declin; - Double_t h0 = hangle*360./24.; //deg - Double_t RaHOffset = fRepDrive->GetRa() - h0; + + TFile rootfile(skyplotfilename, "RECREATE", + "sky plots in some variations"); + fHistSignif.Write(); + fHistNexcess.Write(); + fHistOn.Write(); + fHistSignif.Write(); + + // from Wolfgang: + //-------------------------------------------------------------------- + // Dec-Hour lines + Double_t rho, sinrho, cosrho; + 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); + } + + 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); + + *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] + Double_t gridbinning = fGridBinning; + Double_t gridfinebin = fGridFineBin; + Int_t numgridlines = (Int_t)(4.0/gridbinning); + Double_t aberr = 1.07; + Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ; + + Double_t declin, hangle; // [deg], [h] + MStarCamTrans mstarc(*fGeomCam, *fObservatory); + 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! + + TLatex *pix; + + Char_t tit[100]; + Double_t xtxt; + Double_t ytxt; + + Double_t xlo = -534.0 * mmtodeg; + Double_t xup = 534.0 * mmtodeg; + + Double_t ylo = -534.0 * mmtodeg; + Double_t yup = 534.0 * mmtodeg; + + Double_t xx, yy; + + + // direction for the center of the camera + Double_t dec0 = declin; + 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 ------------------------------------ - - const Int_t Ndec = numgridlines; - Double_t dec[Ndec]; - Double_t ddec = gridbinning; - - Int_t nh = (Int_t)(4.0/gridfinebin); - const Int_t Nh = nh+1; - Double_t h[Nh]; - Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926); - if ( dh > 360.0/((Double_t)(Nh-1)) ) - dh = 360.0/((Double_t)(Nh-1)); - -// start to copy - for (Int_t j=0; jSetOptFit(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 ------------------------------------ + + const Int_t Ndec = numgridlines; + Double_t dec[Ndec]; + Double_t ddec = gridbinning; + + Int_t nh = (Int_t)(4.0/gridfinebin); + const Int_t Nh = nh+1; + Double_t h[Nh]; + Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926); + if ( dh > 360.0/((Double_t)(Nh-1)) ) + dh = 360.0/((Double_t)(Nh-1)); + + // start to copy + for (Int_t j=0; j 90.0) continue; - + for (Int_t k=0; kcd(2); + { + Double_t hh0 = h0 *24.0/360.0; + Double_t hx = h[k]*24.0/360.0; + mstarc.Cel0CelToCam(dec0, hh0, dec[j], hx, + xx, yy); + xx = xx * mmtodeg * aberr; + yy = yy * mmtodeg * aberr; + // xh[k] = xx * mmtodeg * aberr; + // yh[k] = yy * mmtodeg * aberr; + xh[k] = cosg * xx + sing * yy; + yh[k] =-sing * xx + cosg * yy; + // xh[k] = cosrho * xx + sinrho * yy; + // yh[k] =-sinrho * xx + cosrho * yy; + + + //gLog << "dec0, h0 = " << dec0 << ", " << h0 + // << " : dec, h, xh, yh = " << dec[j] << ", " + // << h[k] << "; " << xh[k] << ", " << yh[k] << endl; + } + + // c1->cd(2); TGraph * graph1 = new TGraph(Nh,xh,yh); //graph1->SetLineColor(j+1); @@ -1017,7 +1104,7 @@ graph1->Draw("L"); //delete graph1; -// graphvec.push_back(*graph1); -// graphvec[j].Draw("L"); - + // graphvec.push_back(*graph1); + // graphvec[j].Draw("L"); + sprintf(tit,"Dec = %6.2f", dec[j]); xtxt = xlo + (xup-xlo)*0.1; @@ -1028,63 +1115,62 @@ //pix->Draw(""); //delete pix; - - } -//stop copy - //----- lines for fixed hour angle ------------------------------------ - - Int_t ndec1 = (Int_t)(4.0/gridfinebin); - const Int_t Ndec1 = ndec1; - Double_t dec1[Ndec1]; - Double_t ddec1 = gridfinebin; - - const Int_t Nh1 = numgridlines; - Double_t h1[Nh1]; - Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926); - if ( dh1 > 360.0/((Double_t)Nh1) ) - dh1 = 360.0/((Double_t)Nh1); - -// start copy - for (Int_t j=0; j 360.0/((Double_t)Nh1) ) + dh1 = 360.0/((Double_t)Nh1); + + // start copy + for (Int_t j=0; j 90.0) continue; - - Double_t hh0 = h0 *24.0/360.0; - Double_t hhx = h1[k]*24.0/360.0; - mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx, + { + if (fabs(dec1[j]) > 90.0) continue; + + Double_t hh0 = h0 *24.0/360.0; + Double_t hhx = h1[k]*24.0/360.0; + mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx, xx, yy); -// xd[count] = xx * mmtodeg * aberr; -// yd[count] = yy * mmtodeg * aberr; - - xx = xx * mmtodeg * aberr; - yy = yy * mmtodeg * aberr; - xd[count] = cosg * xx + sing * yy; - yd[count] =-sing * xx + cosg * yy; - - //gLog << "dec0, h0 = " << dec0 << ", " << h0 - // << " : dec1, h1, xd, yd = " << dec1[j] << ", " - // << h1[k] << "; " << xd[count] << ", " << yd[count] << endl; - - count++; - } - -// c1->cd(2); + // xd[count] = xx * mmtodeg * aberr; + // yd[count] = yy * mmtodeg * aberr; + + xx = xx * mmtodeg * aberr; + yy = yy * mmtodeg * aberr; + xd[count] = cosg * xx + sing * yy; + yd[count] =-sing * xx + cosg * yy; + + //gLog << "dec0, h0 = " << dec0 << ", " << h0 + // << " : dec1, h1, xd, yd = " << dec1[j] << ", " + // << h1[k] << "; " << xd[count] << ", " << yd[count] << endl; + + count++; + } + + // c1->cd(2); TGraph * graph1 = new TGraph(count,xd,yd); //graph1->SetLineColor(k+1); @@ -1111,42 +1197,41 @@ //pix->Draw(""); //delete pix; - - } - -// c1->cd(2); - 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(.05); - c1->cd(1); - pix->Draw(""); - c1->cd(2); - pix->Draw(""); - c1->cd(3); - pix->Draw(""); - //delete pix; - - 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); - c2->Divide(1,2); - c2->SetBorderMode(0); - c2->cd(1); - fHistSignif.Draw("colz"); - c2->cd(2); - fHistNexcess.Draw("colz"); - c2->Write(); - - rootfile.Close(); - delete c1; - delete c2; - delete pix; - + + } + + // c1->cd(2); + 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(.05); + c1->cd(1); + pix->Draw(""); + c1->cd(2); + pix->Draw(""); + c1->cd(3); + pix->Draw(""); + //delete pix; + + 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); + c2->Divide(1,2); + c2->SetBorderMode(0); + c2->cd(1); + fHistSignif.Draw("colz"); + c2->cd(2); + fHistNexcess.Draw("colz"); + c2->Write(); + + rootfile.Close(); + delete c1; + delete c2; + delete pix; } @@ -1157,5 +1242,4 @@ "all the alpha plots"); - vector ::iterator alpha_iterator; int index = 0; // index of the TH2F histograms Char_t strtitle[100]; @@ -1164,8 +1248,8 @@ Int_t nrow, ncolumn, non; - - alpha_iterator = fHistAlpha.begin(); - - while( alpha_iterator != fHistAlpha.end() ) { + TH1D *alpha_iterator = NULL; + TOrdCollectionIter Next(fHistAlpha); + + while( (alpha_iterator = (TH1D*)Next())) { nrow = index/fHistOn.GetNbinsX() + 1; // row of the histogram (y) @@ -1182,8 +1266,7 @@ xpos, ypos, signif, nex, non); - (*alpha_iterator).SetName(strname); - (*alpha_iterator).SetTitle(strtitle); - (*alpha_iterator).Write(); - alpha_iterator++; + alpha_iterator->SetName(strname); + alpha_iterator->SetTitle(strtitle); + alpha_iterator->Write(); index++; } @@ -1192,140 +1275,5 @@ } -// -------------------------------------------------------------------------- -// -// Draw the histogram -// -void MSkyPlot::Draw(Option_t *opt) -{ -/* - TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); - pad->SetBorderMode(0); - - AppendPad(""); - - pad->Divide(1, 2, 0, 0.03); - - TObject *catalog = GetCatalog(); - - // Initialize upper part - pad->cd(1); - gPad->SetBorderMode(0); - gPad->Divide(3, 1); -*/ -/* - // PAD #1 - pad->GetPad(1)->cd(1); - gPad->SetBorderMode(0); - fHist.GetZaxis()->SetRangeUser(0,fAlphaONMax); - TH1 *h3 = fHist.Project3D("yx_on"); - fHist.GetZaxis()->SetRange(0,9999); - h3->SetDirectory(NULL); - h3->SetXTitle(fHist.GetXaxis()->GetTitle()); - h3->SetYTitle(fHist.GetYaxis()->GetTitle()); - h3->Draw("colz"); - h3->SetBit(kCanDelete); - catalog->Draw("mirror same"); - - // PAD #2 - pad->GetPad(1)->cd(2); - gPad->SetBorderMode(0); - fHist.GetZaxis()->SetRange(0,0); - TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning.... - fHist.GetZaxis()->SetRange(0,9999); - h4->SetTitle("Significance"); - h4->SetDirectory(NULL); - h4->SetXTitle(fHist.GetXaxis()->GetTitle()); - h4->SetYTitle(fHist.GetYaxis()->GetTitle()); - h4->Draw("colz"); - h4->SetBit(kCanDelete); - catalog->Draw("mirror same"); - - // PAD #3 - pad->GetPad(1)->cd(3); - gPad->SetBorderMode(0); - fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2); - TH1 *h2 = fHist.Project3D("yx_off"); - h2->SetDirectory(NULL); - h2->SetXTitle(fHist.GetXaxis()->GetTitle()); - h2->SetYTitle(fHist.GetYaxis()->GetTitle()); - h2->Draw("colz"); - h2->SetBit(kCanDelete); - catalog->Draw("mirror same"); - - - // Initialize lower part - pad->cd(2); - gPad->SetBorderMode(0); - gPad->Divide(3, 1); - - // PAD #4 - pad->GetPad(2)->cd(1); - gPad->SetBorderMode(0); - TH1 *h1 = fHist.ProjectionZ("Alpha"); - h1->SetDirectory(NULL); - h1->SetTitle("Distribution of \\alpha"); - h1->SetXTitle(fHist.GetZaxis()->GetTitle()); - h1->SetYTitle("Counts"); - h1->Draw(opt); - h1->SetBit(kCanDelete); - - // PAD #5 - pad->GetPad(2)->cd(2); - gPad->SetBorderMode(0); - TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff"); - h5->Add(h2, -1); - h5->SetTitle("Difference of on- and off-distribution"); - h5->SetDirectory(NULL); - h5->Draw("colz"); - h5->SetBit(kCanDelete); - catalog->Draw("mirror same"); - - // PAD #6 - pad->GetPad(2)->cd(3); - gPad->SetBorderMode(0); - TH1 *h0 = fHist.Project3D("yx_all"); - h0->SetDirectory(NULL); - h0->SetXTitle(fHist.GetXaxis()->GetTitle()); - h0->SetYTitle(fHist.GetYaxis()->GetTitle()); - h0->Draw("colz"); - h0->SetBit(kCanDelete); - catalog->Draw("mirror same"); -*/ -} - -// -------------------------------------------------------------------------- -// -// Everything which is in the main pad belongs to this class! -// -Int_t MSkyPlot::DistancetoPrimitive(Int_t px, Int_t py) -{ - return 0; -} - -// -------------------------------------------------------------------------- -// -// Set all sub-pads to Modified() -// -/* -void MSkyPlot::Modified() -{ - if (!gPad) - return; - - TVirtualPad *padsave = gPad; - padsave->Modified(); - padsave->GetPad(1)->cd(1); - gPad->Modified(); - padsave->GetPad(1)->cd(3); - gPad->Modified(); - padsave->GetPad(2)->cd(1); - gPad->Modified(); - padsave->GetPad(2)->cd(2); - gPad->Modified(); - padsave->GetPad(2)->cd(3); - gPad->Modified(); - gPad->cd(); -} -*/ + // -------------------------------------------------------------------------- // Index: trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h (revision 6189) +++ trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h (revision 6190) @@ -6,6 +6,4 @@ #endif -#include - #ifndef ROOT_TH2 #include @@ -14,4 +12,8 @@ #ifndef ROOT_TH1 #include +#endif + +#ifndef ROOT_TOrdCollection +#include #endif @@ -29,24 +31,29 @@ class MHillasSrc; class MNewImagePar; +class MHadronness; class MGeomCam; - +class TOrdCollection; class MSkyPlot : public MTask { private: - MGeomCam *fGeomCam; //! container to take the event time from + + MGeomCam *fGeomCam; //! container to take the event time from MTime *fTime; //! container to take the event time from MPointingPos *fPointPos; //! container to take pointing position from - MReportDrive *fRepDrive; - MSrcPosCam *fSrcPosCam; //! container with x and y of the source - MSrcPosCam *fPntPosCam; //! container with x and y of the position MReportDrive.GetRa, MReportDrive.GetDec + MReportDrive *fRepDrive; //! + MSrcPosCam *fSrcPosCam; //! container with x and y of the source + MSrcPosCam *fPntPosCam; //! container with x and y of the position MReportDrive.GetRa, MReportDrive.GetDec MObservatory *fObservatory; //! container to take observatory location from - MHillas *fHillas; - MHillasExt *fHillasExt; - MHillasSrc *fHillasSrc; - MNewImagePar *fNewImagePar; + MHillas *fHillas; //! + MHillasExt *fHillasExt; //! + MHillasSrc *fHillasSrc; //! + MNewImagePar *fNewImagePar; //! + MHadronness *fHadron; //! - Float_t fMm2Deg; // conversion factor for display in degrees - Double_t fGridBinning; // degrees - Double_t fGridFineBin; // degrees + TOrdCollection *fHistAlpha; // vector of histograms for alpha + + Float_t fMm2Deg; // conversion factor for display in degrees + Double_t fGridBinning; // degrees + Double_t fGridFineBin; // degrees // Float_t fAlphaCut; // Alpha cut @@ -59,5 +66,4 @@ // Float_t fMaxLD; // Maximum distance in percent of dist - std::vector fHistAlpha; // vector of histograms for alpha Int_t fNumalphahist; // number of histograms for alpha Int_t fNumBinsAlpha; @@ -65,28 +71,27 @@ Float_t fAlphaLeftEdge; Float_t fAlphaRightEdge; - Float_t fAlphaONMax; // [deg] , upper cut for alpha ON region in the alpha plot - // [deg], ON region in the alpha plot, maybe 5 deg is better - // NOTE: up to now only values of 5, 10, 15, 20 degrees are possible + Float_t fAlphaONMax; // [deg] , upper cut for alpha ON region in the alpha plot, [deg], ON region in the alpha plot, maybe 5 deg is better, NOTE: up to now only values of 5, 10, 15, 20 degrees are possible Float_t fAlphaBgLow; // lower limit for bg region in the ON alpha plot Float_t fAlphaBgUp; // upper limit for bg region in the ON alpha plot - TH2D fHistSignif; // sky plot of significance vs. x and y - 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; - Double_t fDec0; - Bool_t kSaveAlphaPlots; - Bool_t kSaveSkyPlots; - Bool_t kSaveNexPlot; - - Float_t fMinXGrid; // [deg] , left edge of the skyplot - Float_t fMaxXGrid; // [deg] , right edge of the skyplot - Float_t fMinYGrid; // [deg] , upper edge of the skyplot - Float_t fMaxYGrid; // [deg] , lower edge of the skyplot - Float_t fBinStepGrid; // [deg] , grid size - Int_t fNumStepsX; // number of bins in x - Int_t fNumStepsY; // number of bins in y + TH2D fHistSignif; // sky plot of significance vs. x and y + 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 + Bool_t fUseRF; // use RF hadronness cut instead of supercuts + Double_t fRa0; // + Double_t fDec0; // + Bool_t fSaveAlphaPlots; // + Bool_t fSaveSkyPlots; // + Bool_t fSaveNexPlot; // + + Float_t fMinXGrid; // [deg] , left edge of the skyplot + Float_t fMaxXGrid; // [deg] , right edge of the skyplot + Float_t fMinYGrid; // [deg] , upper edge of the skyplot + Float_t fMaxYGrid; // [deg] , lower edge of the skyplot + Float_t fBinStepGrid; // [deg] , grid size + Int_t fNumStepsX; // number of bins in x + Int_t fNumStepsY; // number of bins in y @@ -97,4 +102,5 @@ Float_t fMaxDist; // dist max cut (ever possible) Float_t fMinDist; // dist min cut (ever possible) + Float_t fHadrCut; // hadronness cut // coefficients for the cuts: @@ -107,6 +113,5 @@ TString fSkyHName; // name for histogram with sky plots TString fNexHName; // name for canvas with Nex plot - - Int_t DistancetoPrimitive(Int_t px, Int_t py); + Float_t fElaTime; // elapsed time [min] TObject *GetCatalog(); @@ -119,10 +124,26 @@ public: MSkyPlot(const char *name=NULL, const char *title=NULL); + ~MSkyPlot(); -// TH1 *GetHistByName(const TString name) { return &fHist; } - TH2D *GetHistSignif() { return &fHistSignif; } - TH2D *GetHistNexcess() { return &fHistNexcess; } - TH2D *GetHistOn() { return &fHistOn; } + Double_t CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2); + + TH2D *GetHistSignif () { return &fHistSignif; } + TH2D *GetHistNexcess () { return &fHistNexcess; } + TH2D *GetHistOn () { return &fHistOn; } TH1D *GetHistSignifGaus() { return &fHistSignifGaus; } + + Int_t GetSize() const { return fHistAlpha->GetSize(); } + + TH1D *GetAlphaPlot( const Int_t i=-1); + const TH1D *GetAlphaPlot( const Int_t=-1) const; + + void ReadCuts(const TString parSCinit); + + void SaveAlphaPlots(const TString stri2); + void SaveNexPlot(const TString stri3); + void SaveSkyPlots(TString stri); + + void SetAlphaCut(Float_t alpha); + void SetAlphaBGLimits(Float_t alphalow, Float_t alphalup); void SetMinDist(Float_t dist) { fMinDist = dist; } // Absolute minimum distance @@ -130,31 +151,20 @@ void SetSizeMin(Float_t size) { fSizeMin = size; } // Absolute minimum Size void SetSizeMax(Float_t size) { fSizeMax = size; } // Absolute maximum Size - void ReadCuts(const TString parSCinit); 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 SetHadrCut(Float_t b) { fHadrCut = b; } // hadronness cut void SetOutputSkyName(TString outname2) { fSkyHName = outname2; } - void SaveSkyPlots(TString stri); - void SetNotSaveSkyPlots() { kSaveSkyPlots = kFALSE; } - + void SetNotSaveSkyPlots() { fSaveSkyPlots = kFALSE; } void SetOutputAlphaName(TString outname1) { fAlphaHName = outname1; } - void SaveAlphaPlots(const TString stri2); - void SetNotSaveAlphaPlots() { kSaveAlphaPlots = kFALSE; } + void SetNotSaveAlphaPlots() { fSaveAlphaPlots = kFALSE; } void SetOutputNexName(TString outname3) { fNexHName = outname3; } - void SaveNexPlot(const TString stri3); - void SetNotSaveNexPlot() { kSaveNexPlot = kFALSE; } + void SetElapsedTime(Float_t g) { fElaTime = g; } + void SetNotSaveNexPlot() { fSaveNexPlot = kFALSE; } + void SetUseRF() { fUseRF = kTRUE; } - void SetAlphaCut(Float_t alpha); - void SetAlphaBGLimits(Float_t alphalow, Float_t alphalup); - - //std::vector ::iterator GetAlphaPlots() { std::vector ::iterator iter; return iter = fHistAlpha.begin() ; } - std::vector GetAlphaPlots() { return fHistAlpha ; } - - void Draw(Option_t *option=""); - - ClassDef(MSkyPlot, 0) //2D-histogram in alpha, x and y + ClassDef(MSkyPlot, 1) //2D-histogram in alpha, x and y };