Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTime.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTime.cc (revision 4409) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTime.cc (revision 4410) @@ -69,5 +69,5 @@ totalLiveTime += fLiveTimeBin[bin]; } - *fLog << GetName() << ": Total live time " << totalLiveTime << endl; + *fLog << GetName() << ": Total live time " << totalLiveTime << " sec" << endl; } else Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPlace.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPlace.h (revision 4409) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPlace.h (revision 4410) @@ -25,5 +25,5 @@ TH2F* fHistPos; // histogram of the used source positions - OnOffMode_t fMode; // On/Off data mode (write/read to/from the histogram) + OnOffMode_t fMode; // On/Off data mode (write/read to/from the histogram) Bool_t fCreateHisto; // flag to decide whether internal histo is created or not @@ -48,4 +48,5 @@ void SetInternalHistoName(TString name) {fHistoName=name;} void SetInternalHistoBinSize(Float_t size){fHistoBinPrec=size;} + void SetPositionHisto(TH2F* histo) {fHistPos=histo;} void SetCreateHisto(Bool_t inp=kTRUE) {fCreateHisto=inp;} Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.cc (revision 4409) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.cc (revision 4410) @@ -40,4 +40,9 @@ #include +#include +#include +#include +#include +#include #include "MSrcPosCam.h" @@ -61,9 +66,11 @@ // Default constructor. // -MSrcPosFromStars::MSrcPosFromStars(Float_t first, Float_t second, const char *name, const char *title) - : fStars(NULL), fNumStars(2), fDistanceFirstStar(first), fDistanceSecondStar(second) +MSrcPosFromStars::MSrcPosFromStars(const char *name, const char *title) + : fStars(NULL) { fName = name ? name : gsDefName.Data(); fTitle = title ? title : gsDefTitle.Data(); + + fDistances.Set(0); } @@ -84,4 +91,6 @@ } + fNumStars = fDistances.GetSize(); + return kTRUE; } @@ -90,8 +99,32 @@ // // +static Bool_t HandleInput() +{ + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + while (1) + { + // + // While reading the input process gui events asynchronously + // + timer.TurnOn(); + cout << "Type 'q' to exit, to go on: " << endl; + TString input; + input = getchar(); + timer.TurnOff(); + + if (input=="q\n") + return kFALSE; + + if (input=="\n") + return kTRUE; + }; + + return kFALSE; +} + Int_t MSrcPosFromStars::ComputeNewSrcPosition() { - if (fDistanceFirstStar == 0. || fDistanceSecondStar == 0.) + if (fNumStars == 0) { if (fStars->GetNumStars() > 0) @@ -135,50 +168,167 @@ } } - else if (fStars->GetNumStars() == fNumStars) + else if (fStars->GetNumStars() >= fNumStars) { - MStarLocalPos& firstStar = (*fStars)[0]; - MStarLocalPos& secondStar = (*fStars)[1]; - - if (firstStar.GetChiSquareNdof() > 0. && firstStar.GetChiSquareNdof() < 10. && - secondStar.GetChiSquareNdof() > 0. && secondStar.GetChiSquareNdof() < 10.) + Float_t diameterInnerPixel = 30; //[mm] + Float_t diameterOuterPixel = 60; //[mm] + + Double_t probability = 1.; + + // Resolution and computing time are proportional to bins^2 + const Int_t bins = 200; + Double_t max_x_mm = 600; + Double_t min_x_mm = -max_x_mm; + Double_t max_y_mm = max_x_mm; + Double_t min_y_mm = -max_x_mm; + + // bins to mmrees + const Double_t bin2mm = 2 * max_x_mm / bins; + + // First run, to find the maximum peak and define the area + TH2F *hgrid = new TH2F("hgrid", "", bins, min_x_mm, max_x_mm, bins, min_y_mm, max_y_mm); + + for (Int_t ny = 1; ny <= bins; ny++) + for (Int_t nx = 1; nx <= bins; nx++) + hgrid->SetBinContent(nx, ny, 1.); + + for (UInt_t numstar = 0; numstar < fNumStars; numstar++) { - - Float_t firstStarPosX = firstStar.GetMeanX(); - Float_t firstStarPosY = firstStar.GetMeanY(); - - Float_t secondStarPosX = secondStar.GetMeanX(); - Float_t secondStarPosY = secondStar.GetMeanY(); - - Float_t distanceStars = sqrt(pow(firstStarPosX - secondStarPosX,2) + pow(firstStarPosY - secondStarPosY,2)); - Float_t sin_alpha = (secondStarPosY - firstStarPosY) / distanceStars; - Float_t cos_alpha = (secondStarPosX - firstStarPosX) / distanceStars; - - Float_t x = (pow(fDistanceFirstStar,2) - pow(fDistanceSecondStar,2) + pow(distanceStars,2)) / (2 * distanceStars); - - Float_t arg = 4 * pow(distanceStars,2) * pow(fDistanceFirstStar,2) - pow(pow(fDistanceFirstStar,2) - pow(fDistanceSecondStar,2) + pow(distanceStars,2),2); - - if (arg >= 0.) - { - - Float_t y = sqrt(arg) / (2 * distanceStars); - - Float_t xc1 = firstStarPosX + x * cos_alpha - y * sin_alpha; - Float_t yc1 = firstStarPosY + x * sin_alpha + y * cos_alpha; - - Float_t xc2 = firstStarPosX + x * cos_alpha + y * sin_alpha; - Float_t yc2 = firstStarPosY + x * sin_alpha - y * cos_alpha; - - if (sqrt(xc1*xc1 + yc1*yc1) < sqrt(xc2*xc2 + yc2*yc2)) - GetOutputSrcPosCam()->SetXY(xc1,yc1); - else - GetOutputSrcPosCam()->SetXY(xc2,yc2); - } - else - *fLog << warn << GetName() << " negative argument [" << arg << "] for sqrt()" << endl; - + probability = 1; + + MStarLocalPos& star = (*fStars)[numstar]; + + if (star.GetChiSquareNdof() > 0. && star.GetChiSquareNdof() < 10.) + { + + Float_t starPosX = star.GetMeanX(); + Float_t starPosY = star.GetMeanY(); + Float_t starDist = Dist(0.,0.,starPosX,starPosY); + Float_t starSigma = (starDist<350.?diameterInnerPixel:diameterOuterPixel); + +// *fLog << dbg << "Star[" << numstar << "] pos (" << starPosX << "," << starPosY << ") dist " << starDist << " size " << starSigma << endl; + + if (starSigma != 0) + { + for (Int_t ny = 1; ny < bins + 1; ny++) + { + for (Int_t nx = 0; nx < bins + 1; nx++) + { + Float_t dist = Dist(min_x_mm + nx * bin2mm, starPosX, min_y_mm + ny * bin2mm, starPosY); + Float_t prob = Prob(dist, fDistances[numstar], starSigma); + +// if ( prob > 0.8) +// *fLog << dbg << " x " << min_x_mm + nx * bin2mm << " y " << min_y_mm + ny * bin2mm << " dist " << dist << " stardist " << fDistances[numstar] << " prob " << prob << endl; + + probability = hgrid->GetBinContent(nx, ny)*prob; + hgrid->SetBinContent(nx, ny, probability); + + } + } + } + else probability *= 1; + + } } + + // Finding the peak + Double_t peak[2] = {0,0}; + Double_t peak_value = 0; + + for (Int_t ny = 0; ny < bins + 1; ny++) + { + for (Int_t nx = 0; nx < bins + 1; nx++) + { + if (hgrid->GetBinContent(nx, ny) > peak_value) + { + peak_value = hgrid->GetBinContent(nx, ny); + peak[0] = min_x_mm + nx * bin2mm; + peak[1] = min_y_mm + ny * bin2mm; + } + } + } + + *fLog << dbg << "The peak is at (x, y) = (" << peak[0] << ", " << peak[1] << ") mm" << endl; + + +// TCanvas c1; +// hgrid->Draw("lego"); +// if(!HandleInput()) +// exit(1); + + + // Here we define a small area surrounding the peak to avoid wasting time and resolution anywhere else + + min_x_mm = peak[0] - diameterInnerPixel; + max_x_mm = peak[0] + diameterInnerPixel; + min_y_mm = peak[1] - diameterInnerPixel; + max_y_mm = peak[1] + diameterInnerPixel; + + const Double_t xbin2mm = (max_x_mm - min_x_mm) / bins; + const Double_t ybin2mm = (max_y_mm - min_y_mm) / bins; + + // Other run centered in the peak for more precision + TH2F *hagrid = new TH2F("hagrid", "", bins, min_x_mm, max_x_mm, bins, min_y_mm, max_y_mm); + + for (Int_t ny = 1; ny <= bins; ny++) + for (Int_t nx = 1; nx <= bins; nx++) + hagrid->SetBinContent(nx, ny, 1.); + + + for (UInt_t numstar = 0; numstar < fNumStars; numstar++) + { + probability = 1; + + MStarLocalPos& star = (*fStars)[numstar]; + + if (star.GetChiSquareNdof() > 0. && star.GetChiSquareNdof() < 10.) + { + + Float_t starPosX = star.GetMeanX(); + Float_t starPosY = star.GetMeanY(); + Float_t starDist = Dist(0.,0.,starPosX,starPosY); + Float_t starSigma = (starDist<350.?diameterInnerPixel:diameterOuterPixel); + + if (starSigma != 0) + { + for (Int_t ny = 0; ny < bins; ny++) + { + for (Int_t nx = 0; nx < bins; nx++) + { + Float_t prob = Prob(Dist(min_x_mm + nx * xbin2mm, starPosX, min_y_mm + ny * ybin2mm, starPosY), fDistances[numstar], starSigma); + + probability = hagrid->GetBinContent(nx, ny)*prob; + hagrid->SetBinContent(nx, ny, probability); + } + } + } + else probability *= 1; + + } + } + + // This time we fit the histogram with a 2D gaussian + // Although we don't expect our function to be gaussian... + TF2 *gauss2d = new TF2("gauss2d","[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))*exp(-(y-[3])*(y-[3])/(2*[4]*[4]))", min_x_mm, max_x_mm, min_y_mm, max_y_mm); + + gauss2d->SetParName(0,"f0"); + gauss2d->SetParName(1,"meanx"); + gauss2d->SetParName(2,"sigmax"); + gauss2d->SetParName(3,"meany"); + gauss2d->SetParName(4,"sigmay"); + + gauss2d->SetParameter(0,10); + gauss2d->SetParameter(1,peak[0]); + gauss2d->SetParameter(2,0.002); + gauss2d->SetParameter(3,peak[1]); + gauss2d->SetParameter(4,0.002); + + GetOutputSrcPosCam()->SetXY(gauss2d->GetParameter(1), gauss2d->GetParameter(3)); + + delete hgrid; + delete hagrid; } + return kTRUE; } Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.h (revision 4409) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.h (revision 4410) @@ -1,4 +1,12 @@ #ifndef MARS_MSrcPosFromStars #define MARS_MSrcPosFromStars + +#ifndef ROOT_TArrayF +#include "TArrayF.h" +#endif + +#ifndef ROOT_TMath +#include "TMath.h" +#endif #ifndef MARS_MSrcPlace @@ -15,7 +23,6 @@ MStarLocalCam *fStars; - const UInt_t fNumStars; - Float_t fDistanceFirstStar; - Float_t fDistanceSecondStar; + UInt_t fNumStars; + TArrayF fDistances; virtual Int_t ComputeNewSrcPosition(); @@ -24,8 +31,19 @@ public: - MSrcPosFromStars(Float_t first=0., Float_t second=0., const char *name=NULL, const char *title=NULL); + MSrcPosFromStars(const char *name=NULL, const char *title=NULL); + + +// This fuction is to compute the probability using a ring distribution + Double_t Prob(Double_t d, Double_t dist, Double_t sigma) + { + return TMath::Exp(-(d-dist)*(d-dist)/(2*sigma*sigma));///(sigma*TMath::Sqrt(2*3.141592654)); + } + + Double_t Dist(Double_t x1, Double_t x2, Double_t y1, Double_t y2) + { + return TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); + } - void SetDistanceFirstStar(Float_t dist) { fDistanceFirstStar = dist; } - void SetDistanceSecondStar(Float_t dist) { fDistanceSecondStar = dist; } + void SetDistances(TArrayF dist) { fDistances = dist; } ClassDef(MSrcPosFromStars, 0) // task to calculate the position of the source using the position of stars Index: trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile (revision 4409) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile (revision 4410) @@ -66,5 +66,5 @@ MDispCalc.cc \ MControlPlots.cc \ - MSrcPosFromStars.c \ + MSrcPosFromStars.cc \ MLiveTime.cc \ MLiveTimeCalc.cc