Changeset 7365 for trunk/MagicSoft/Mars/mmuon

Timestamp:

09/26/05 19:17:18 (19 years ago)

Author:

meyer

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mmuon

Files:

• MHMuonPar.cc (modified) (2 diffs)
• MHSingleMuon.cc (modified) (9 diffs)
• MHSingleMuon.h (modified) (2 diffs)
• MMuonSearchPar.cc (modified) (3 diffs)
• MMuonSearchPar.h (modified) (2 diffs)
• MMuonSetup.cc (modified) (1 diff)

trunk/MagicSoft/Mars/mmuon/MHMuonPar.cc

@@ -258 +258 @@
     const Int_t bin2 = p.GetXaxis()->FindFixBin(b);
 
-    return Integral(p, bin1, bin2);
+//    return Integral(p, bin1, bin2);
+    return p.Integral(bin1, bin2);
 }
 
@@ -265 +266 @@
     if (TString(opt)==TString("pad4"))
     {
-        const TString txt = Form("\\Sigma_{%.2f\\circ}^{%.2f\\circ} = %3.1fm",
-                                 fgIntegralLoLim, fgIntegralUpLim, Integral(fHistBroad)*1000);
+        const TString txt = Form("\\Sigma_{%.2f\\circ}^{%.2f\\circ} = %.3f",
+                                 fgIntegralLoLim, fgIntegralUpLim, Integral(fHistBroad));
+//                                 fgIntegralLoLim, fgIntegralUpLim, Integral(fHistBroad)*1000);
 
         TLatex text(0.55, 0.93, txt);

trunk/MagicSoft/Mars/mmuon/MHSingleMuon.cc

@@ -127 +127 @@
     fHistWidth.UseCurrentStyle();
 
+    fHistTime.SetName("HistTime");
+    fHistTime.SetTitle("HistTime");
+    fHistTime.SetXTitle("timing difference");
+    fHistTime.SetYTitle("number of pixels");
+    fHistTime.SetDirectory(NULL);
+    fHistTime.SetFillStyle(4000);
+    fHistTime.UseCurrentStyle();
+
     MBinning bins;
     bins.SetEdges(20, -180, 180);
@@ -133 +141 @@
     bins.SetEdges(28, 0.3, 1.7);
     bins.Apply(fHistWidth);
+
+    bins.SetEdges(100, -10, 10);
+    bins.Apply(fHistTime);
 }
 
@@ -169 +180 @@
     ApplyBinning(*plist, "ArcPhi",    &fHistPhi);
     ApplyBinning(*plist, "MuonWidth", &fHistWidth);
+    ApplyBinning(*plist, "MuonTime", &fHistTime);
 
     return kTRUE;
@@ -182 +194 @@
     fHistPhi.Reset();
     fHistWidth.Reset();
+    fHistTime.Reset();
 
     const Int_t entries = fSignalCam->GetNumPixels();
@@ -202 +215 @@
         if (dist < fMuonSearchPar->GetRadius() + fMargin &&
             dist > fMuonSearchPar->GetRadius() - fMargin)
-            fHistPhi.Fill(TMath::ATan2(dx, dy)*TMath::RadToDeg(), pix.GetNumPhotons());
-
+        {
+            fHistTime.Fill(pix.GetArrivalTime()-fMuonSearchPar->GetTime());
+        }
         // use only the inner pixles. FIXME: This is geometry dependent
         if(i>397)
@@ -210 +224 @@
         fHistWidth.Fill(dist*fGeomCam->GetConvMm2Deg(), pix.GetNumPhotons());
     }
-
+    // Setup the function and perform the fit
+    TF1 g1("g1", "gaus", -10, 10);
+
+    // Choose starting values as accurate as possible
+    g1.SetParameter(0, fHistTime.GetMaximum());
+    g1.SetParameter(1, 0);
+    g1.SetParameter(2, 0.3);
+    g1.SetParLimits(1, -0.5, 0.5);
+    g1.SetParLimits(2, 0, 1);
+    // options : N  do not store the function, do not draw
+    //           I  use integral of function in bin rather than value at bin center
+    //           R  use the range specified in the function range
+    //           Q  quiet mode
+    fHistTime.Fit(&g1, "QNRB");
+
+    //    Double_t err;
+    Double_t mean, meanerr, sig, sigerr;
+    gMinuit->GetParameter(2, sig, sigerr); // get the sigma value
+    gMinuit->GetParameter(1, mean, meanerr); // get the sigma value
+
+    for (Int_t i=0; i<entries; i++)
+    {
+        const MSignalPix &pix  = (*fSignalCam)[i];
+        const MGeomPix   &gpix = (*fGeomCam)[i];
+
+        const Float_t dx = gpix.GetX() - cenx;
+        const Float_t dy = gpix.GetY() - ceny;
+
+        const Float_t dist = TMath::Hypot(dx, dy);
+
+        // if the signal is not near the estimated circle, it is ignored.
+        if (dist < fMuonSearchPar->GetRadius() + fMargin &&
+            dist > fMuonSearchPar->GetRadius() - fMargin)
+        {
+            if(TMath::Abs(pix.GetArrivalTime()-(fMuonSearchPar->GetTime()+mean))<2*sig)
+                fHistPhi.Fill(TMath::ATan2(dx, dy)*TMath::RadToDeg(), pix.GetNumPhotons());
+        }
+    }
     // error estimation (temporaly)
     //  The error is estimated from the signal. In order to do so, we have to
@@ -244 +295 @@
 // below the threshold is found, kFALSE is returned.
 //
-Bool_t MHSingleMuon::FindRangeAboveThreshold(const TH1 &h, Float_t thres, Int_t &first, Int_t &last) const
+Bool_t MHSingleMuon::FindRangeAboveThreshold(const TProfile &h, Float_t thres, Int_t &first, Int_t &last) const
 {
     const Int_t n      = h.GetNbinsX();
@@ -362 +413 @@
     //           R  use the range specified in the function range
     //           Q  quiet mode
-    fHistWidth.Fit(&f1, "QR0");
+//    fHistWidth.Fit(&f1, "QRO");
+    fHistWidth.Fit(&f1, "QRN");
 
     chi = f1.GetChisquare()/f1.GetNDF();
@@ -477 +529 @@
     gPad->SetBorderMode(0);
     fHistWidth.Draw();
-
-}
+}

trunk/MagicSoft/Mars/mmuon/MHSingleMuon.h

@@ -24 +24 @@
     Double_t fMargin;               //!
 
-    TH1F     fHistPhi;   // Histogram of photon distribution along the arc.
+    TProfile fHistPhi;   // Histogram of photon distribution along the arc.
     TProfile fHistWidth; // Histogram of radial photon distribution of the arc.
+    TH1F     fHistTime;   // Histogram of arrival time distribution along the arc.
 
-    Bool_t FindRangeAboveThreshold(const TH1 &h, Float_t thres, Int_t &first, Int_t &last) const;
+    Bool_t FindRangeAboveThreshold(const TProfile &h, Float_t thres, Int_t &first, Int_t &last) const;
 
 public:
@@ -38 +39 @@
     Bool_t CalcWidth(Double_t, Double_t &, Double_t &);
 
-    const TH1F &GetHistPhi() const { return fHistPhi; }
+    const TProfile &GetHistPhi() const { return fHistPhi; }
+    const TH1F     &GetHistTime() const { return fHistTime; }
     const TProfile &GetHistWidth() const { return fHistWidth; }
 

trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc

@@ -89 +89 @@
     fCenterX   =  0;
     fCenterY   =  0;
+    fTime      =  0;
 }
 
@@ -150 +151 @@
 
     Int_t p=0;
+    Int_t q=0;
+
+    Float_t fTime2=0;
+
     for (int i=0; i<n; i++)
     {
         const MSignalPix &pix = evt[i];
+
         if (pix.IsPixelUsed())
         {
@@ -160 +166 @@
             fY[p] = geom[i].GetY();
             p++;
+
+            //timing
+            if(pix.IsPixelCore())
+            {
+                fTime  += pix.GetArrivalTime();
+                fTime2 += pix.GetArrivalTime()*pix.GetArrivalTime();
+                q++;
+            }
         }
     }
+    fTime    = fTime/q;
+    fTime2   = fTime2/q;
+    fTimeRms = TMath::Sqrt(fTime2-(fTime*fTime));
+
     fSignal.Set(p);
 

trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h

@@ -21 +21 @@
     Float_t fCenterX;   // An estimated center position in X of the muon ring [mm]
     Float_t fCenterY;   // An estimated center position in Y of the muon ring [mm]
+    Float_t fTime;
+    Float_t fTimeRms;
 
     MArrayF fSignal;    //! Temporary storage for signal
@@ -40 +42 @@
     Float_t GetCenterX()   const { return fCenterX; }
     Float_t GetCenterY()   const { return fCenterY; }
+    Float_t GetTime()      const { return fTime; }
+    Float_t GetTimeRms()   const { return fTimeRms; }
 
     // MMuonSearchPar

trunk/MagicSoft/Mars/mmuon/MMuonSetup.cc

@@ -49 +49 @@
 //
 MMuonSetup::MMuonSetup(const char *name, const char *title)
-    : fMargin(0.2), fThresholdArcPhi(30), fThresholdArcWidth(2)
+    : fMargin(0.2), fThresholdArcPhi(5), fThresholdArcWidth(2)
 {
     fName  = name  ? name  : "MMuonSetup";