Changeset 4410 for trunk/MagicSoft/Mars

Timestamp:

07/20/04 12:15:06 (20 years ago)

Author:

jlopez

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mtemp/mifae/library

Files:

• MLiveTime.cc (modified) (1 diff)
• MSrcPlace.h (modified) (2 diffs)
• MSrcPosFromStars.cc (modified) (5 diffs)
• MSrcPosFromStars.h (modified) (3 diffs)
• Makefile (modified) (1 diff)

trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTime.cc

@@ -69 +69 @@
             totalLiveTime += fLiveTimeBin[bin];
         }               
-        *fLog << GetName() << ": Total live time " << totalLiveTime << endl;
+        *fLog << GetName() << ": Total live time " << totalLiveTime << " sec" << endl;
     }
     else

trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPlace.h

@@ -25 +25 @@
 
   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 +48 @@
   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;}
 

trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.cc

@@ -40 +40 @@
 
 #include <TList.h>
+#include <TH2F.h>
+#include <TF2.h>
+#include <TTimer.h>
+#include <TString.h>
+#include <TCanvas.h>
 
 #include "MSrcPosCam.h"
@@ -61 +66 @@
 // 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 +91 @@
     }
   
+  fNumStars = fDistances.GetSize();
+  
   return kTRUE;
 }
@@ -90 +99 @@
 //
 // 
+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, <return> 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 +168 @@
         }
     }
-  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;
 }

trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.h

@@ -1 +1 @@
 #ifndef MARS_MSrcPosFromStars
 #define MARS_MSrcPosFromStars
+
+#ifndef ROOT_TArrayF
+#include "TArrayF.h"
+#endif
+
+#ifndef ROOT_TMath
+#include "TMath.h"
+#endif
 
 #ifndef MARS_MSrcPlace
@@ -15 +23 @@
   MStarLocalCam *fStars;
 
-  const UInt_t fNumStars;
-  Float_t fDistanceFirstStar;
-  Float_t fDistanceSecondStar;
+  UInt_t fNumStars;
+  TArrayF fDistances;
 
   virtual Int_t ComputeNewSrcPosition();
@@ -24 +31 @@
 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

trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile

@@ -66 +66 @@
         MDispCalc.cc \
         MControlPlots.cc \
-        MSrcPosFromStars.c \
+        MSrcPosFromStars.cc \
         MLiveTime.cc \
         MLiveTimeCalc.cc