Changeset 1888 for trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc

Timestamp:

04/02/03 09:03:22 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (modified) (10 diffs)

trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc

@@ -79 +79 @@
 #include "MMcEvt.hxx"
 #include "MMcTrig.hxx"
+#include "MBinning.h"
+
+#include "TRandom3.h"
+#include "MParameters.h"
 
 ClassImp(MCT1ReadPreProc);
@@ -119 +123 @@
 // Add this file as the last entry in the chain
 //
-Int_t MCT1ReadPreProc::AddFile(const char *txt, Int_t)
+void MCT1ReadPreProc::AddFile(const char *txt)
 {
     const char *name = gSystem->ExpandPathName(txt);
@@ -129 +133 @@
     {
         *fLog << warn << "WARNING - Problem reading header... ignored." << endl;
-        return 0;
+        return;
     }
 
@@ -136 +140 @@
     {
         *fLog << warn << "WARNING - File contains no data... ignored." << endl;
-        return 0;
+        return;
     }
 
@@ -144 +148 @@
 
     fFileNames->AddLast(new TNamed(txt, ""));
-    return 1;
 }
 
@@ -703 +706 @@
 
     //
+    //  look for the HourAngle container in the plist
+    //
+    fHourAngle = (MParameterD*)pList->FindCreateObj("MParameterD","HourAngle");
+    if (!fHourAngle)
+        return kFALSE;
+    fHourAngle->SetTitle("Store the CT1 hour angle [deg]");
+
+    //
+    //  look for the ThetaOrig container in the plist
+    //
+    fThetaOrig = (MParameterD*)pList->FindCreateObj("MParameterD","ThetaOrig");
+    if (!fThetaOrig)
+        return kFALSE;
+    fThetaOrig->SetTitle("Store the original CT1 zenith angle [rad]");
+
+    //
     //  look for the MCerPhotEvt class in the plist
     //
@@ -772 +791 @@
 
     return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Smear Theta uniformly in a bin of Theta; result is stored in ThetaSmeared
+// 
+//
+Bool_t MCT1ReadPreProc::SmearTheta(MParList *plist, Float_t *Theta,
+                             Float_t *ThetaSmeared)
+{
+  // both Theta and ThetaSmeared are in [radians]
+  // the edges are in [degrees]
+
+  const MBinning *binstheta = (MBinning*)plist->FindObject("BinningTheta");
+  if (!binstheta)
+  {
+    *fLog << err << dbginf << "BinningTheta not found ... aborting." << endl;
+    return kFALSE;
+  }
+
+  Int_t nedges = binstheta->GetNumEdges();
+  Double_t *edges  = (Double_t*)binstheta->GetEdges();
+
+  Int_t bin = -1;
+  *ThetaSmeared = *Theta;
+
+  Float_t Thetadeg = (*Theta) * 180.0/TMath::Pi(); 
+  Float_t ThetaSmeareddeg;
+ 
+  // search Theta bin
+  Int_t i;
+  for (i=1; i<nedges; i++)
+  {
+    if (Thetadeg >= *(edges+i)  ) continue;
+    if (Thetadeg  < *(edges+i-1)) break;
+    bin = i;
+    break;
+  }
+
+  Float_t low=0.0;
+  Float_t up =0.0;
+
+  // smear Theta within the Theta bin
+  ThetaSmeareddeg = -1.0;
+  if (bin != -1)
+  {
+    low = *(edges+bin-1);
+    up  = *(edges+bin);
+
+    Double_t ran = ran3.Rndm(1);
+    ThetaSmeareddeg = (low + ran * (up-low));
+  }
+  *ThetaSmeared = ThetaSmeareddeg * TMath::Pi()/180.0;
+
+  //*fLog << "SmearTheta : Thetadeg, ThetaSmeareddeg, low, up, bin = " 
+  //      << Thetadeg 
+  //      << ",  " << ThetaSmeareddeg << ",  " << low << ",  " << up << ",  "
+  //      << bin << endl;    
+
+  return kTRUE;
 }
 
@@ -860 +940 @@
     // int ipreproc_alt_arcs; // "should be" alt according to preproc (arcseconds)
     // int ipreproc_az_arcs;  // "should be" az according to preproc (arcseconds)
+
+    // smear Theta in its Theta bin
+    Float_t ThetaOrig =  TMath::Pi()*(0.5-1./180*event.ialt_arcs/3600); // [radians] 
+    Float_t ThetaSmeared;                                               // [radians]
+    SmearTheta(fParList, &ThetaOrig, &ThetaSmeared);
+    fThetaOrig->SetVal(ThetaOrig);
+
+    // store hour angle
+    fHourAngle->SetVal(event.fhourangle);
+
+    //*fLog << "MCt1ReadPreProc::ProcessEvent; fhourangle = " 
+    //      << event.fhourangle << endl;
 
     fMcEvt->Fill(event.isecs_since_midday,     //0, /*fEvtNum*/
@@ -877 +969 @@
                  fIsMcFile ? event.imcimpact_m*100 : 0,
                  TMath::Pi()/180*event.iaz_arcs/3600, // azimuth (arcseconds)
-                 TMath::Pi()*(0.5-1./180*event.ialt_arcs/3600), // altitude (arcseconds)
+                 ThetaSmeared,
                  0, /* fTFirst */
                  0, /* fTLast */
@@ -895 +987 @@
                  0, /* elec */
                  0, /* muon */
-                 event.fhourangle  /* other */
+                 0  /* other */
                  );