Changeset 1888

Timestamp:

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

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• manalysis/MPointingCorr.cc (modified) (3 diffs)
• manalysis/MPointingCorr.h (modified) (2 diffs)
• manalysis/MSelBasic.cc (modified) (8 diffs)
• manalysis/MSelBasic.h (modified) (2 diffs)
• manalysis/MSelFinal.cc (modified) (6 diffs)
• manalysis/MSelFinal.h (modified) (2 diffs)
• manalysis/MSelStandard.cc (modified) (5 diffs)
• manalysis/MSelStandard.h (modified) (2 diffs)
• mfileio/MCT1ReadPreProc.cc (modified) (10 diffs)
• mfileio/MCT1ReadPreProc.h (modified) (6 diffs)
• mfileio/Makefile (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+
+ 2003/04/02: Wolfgang Wittek
+
+   * mfileio/Makefile
+     - mhist added, because MBinning is used in MCT1ReadPreproc
+
+   * mfileio/MCT1ReadPreProc.[h,cc]
+     - new member function SmearTheta
+     - store smeared  Theta in MMcEvt::fTelescopeTheta
+       store original Theta in MParameterD container "ThetaOrig"
+       store fhourangle     in MParameterD container "HourAngle"
+
+   * manalysis/MPointingCorr.[h,cc]
+     - get hour angle from ParameterD container "HourAngle"
+
+   * manalysis/MSelBasic.[h,cc]
+               MSelStandard.[h,cc]
+               MSelFinal.[h,cc]
+     - new member functions SetCuts()
+
 
  2003/04/01: Abelardo Moralejo
@@ -13 +33 @@
 
  2003/03/31: Thomas Bretz
-
-   * mhist/MHArray.[h,cc]:
-     - added default constructor
-     - added Set-function
-     - added Init function
-     - moved code from constructors to Set and Init
  
    * Makefile.conf.linux:

trunk/MagicSoft/Mars/manalysis/MPointingCorr.cc

@@ -38 +38 @@
 #include "MSrcPosCam.h"
 #include "MGeomCam.h"
+#include "MParameters.h"
 
 #include "MLog.h"
@@ -81 +82 @@
     }
 
+    fHourAngle = (MParameterD*)pList->FindObject("HourAngle", "MParameterD");
+    if (!fHourAngle)
+    {
+        *fLog << dbginf << "HourAngle not found... aborting." << endl;
+        return kFALSE;
+    }
+
 
     fSrcPos = (MSrcPosCam*)pList->FindObject(fSrcName, "MSrcPosCam");
@@ -104 +112 @@
    // (cx, cy) is the source position in the camera [mm]
    //
-   Float_t fhourangle = fMcEvt->GetOtherCphFraction();
+   Float_t fhourangle = fHourAngle->GetVal();
+
+   //*fLog << "MPointingCorr::Process; fhourangle = " 
+   //      << fhourangle << endl;
+
    Float_t cx = -0.05132 - 0.001064 * fhourangle 
                          - 3.530e-6 * fhourangle * fhourangle;

trunk/MagicSoft/Mars/manalysis/MPointingCorr.h

@@ -16 +16 @@
 class MMcEvt;
 class MSrcPosCam;
+class MParameterD;
+
 
 class MPointingCorr : public MTask
@@ -23 +25 @@
     MSrcPosCam   *fSrcPos;
     TString       fSrcName;
+    MParameterD  *fHourAngle;
 
     Float_t      fMm2Deg;

trunk/MagicSoft/Mars/manalysis/MSelBasic.cc

@@ -61 +61 @@
     fTitle = title ? title : "Task to evaluate basic cuts";
 
-    ThetaMin =  0.0;
-    ThetaMax = 60.0;
+    fThetaMin =  0.0;
+    fThetaMax = 60.0;
 }
 
@@ -109 +109 @@
     }
 
-    memset(fErrors, 0, sizeof(fErrors));
+    memset(fCut, 0, sizeof(fCut));
 
     return kTRUE;
@@ -123 +123 @@
 Bool_t MSelBasic::Process()
 {
-  /*
-  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-  *fLog << "=========================================================" << endl;
-  *fLog << "" << endl;
-  *fLog << "MmcEvt data : " << endl;
-  *fLog << "" << endl;
-  fMcEvt->Print();
-  *fLog << "" << endl;
-  *fLog << "PartId() = " << fMcEvt->GetPartId() << endl;
-  *fLog << "Energy() = " << fMcEvt->GetEnergy() << endl;
-  *fLog << "Theta()  = " << fMcEvt->GetTheta() << endl;
-
-  *fLog << "Phi()    = " << fMcEvt->GetPhi() << endl;
-  //*fLog << "CoreD()  = " << fMcEvt->GetCoreD() << endl;
-  //*fLog << "CoreX()  = " << fMcEvt->GetCoreX() << endl;
-
-  //*fLog << "CoreY()  = " << fMcEvt->GetCoreY() << endl;
-  *fLog << "Impact() = " << fMcEvt->GetImpact() << endl;
-  //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl;
-
-  //*fLog << "PassPhotAtm() = " << fMcEvt->GetPassPhotAtm() << endl;
-  //*fLog << "PassPhotRef() = " << fMcEvt->GetPassPhotRef() << endl;
-  //*fLog << "PassPhotCone()  = " << fMcEvt->GetPassPhotCone() << endl;
-
-  //*fLog << "PhotElfromShower()    = " << fMcEvt->GetPhotElfromShower() << endl;
-  //*fLog << "PhotElinCamera()  = " << fMcEvt->GetPhotElinCamera() << endl;
-  *fLog << "TelescopePhi()  = " << fMcEvt->GetTelescopePhi() << endl;
-
-  *fLog << "TelescopeTheta()  = " << fMcEvt->GetTelescopeTheta() << endl;
-  *fLog << "Impact() = " << fMcEvt->GetImpact() << endl;
-  //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl;
-  *fLog << "" << endl;
-  *fLog << "=========================================================" << endl;
-
-  *fLog << "=========================================================" << endl;
-  *fLog << "" << endl;
-  *fLog << "MPedestalPix data : " << endl;
-  *fLog << "" << endl;
-
-  Int_t ntot;
-  ntot = fPed->GetSize();
-    *fLog << "MeanRms() :" << endl;    
-  for (Int_t i=0; i<ntot; i++)
-  {
-    MPedestalPix &pix = (*fPed)[i];
-    //*fLog << "Mean()    = " << i << ",  " << pix.GetMean() << endl;    
-    //*fLog << "Sigma()   = " << i << ",  " << pix.GetSigma() << endl;    
-    *fLog << pix.GetMeanRms() << " ";    
-    //*fLog << "SigmaRms()= " << i << ",  " << pix.GetSigmaRms() << endl;    
-  }
-  *fLog << "" << endl;
-
-  *fLog << "" << endl;
-  ntot = fEvt->GetNumPixels();
-    *fLog << "MCerPhotPix :  pix.GetNumPhotons()" << endl;    
-  for (Int_t i=0; i<ntot; i++)
-  {
-    MCerPhotPix &pix = (*fEvt)[i];
-    *fLog << pix.GetNumPhotons() << " ";    
-  }
-  *fLog << "" << endl;
-
-  *fLog << "" << endl;
-  ntot = fEvt->GetNumPixels();
-    *fLog << "MCerPhotPix :  pix.GetErrorPhot()" << endl;    
-  for (Int_t i=0; i<ntot; i++)
-  {
-    MCerPhotPix &pix = (*fEvt)[i];
-    *fLog << pix.GetErrorPhot() << " ";    
-  }
-  *fLog << "" << endl;
-
-  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-  */
-
-
     Int_t rc = 0;
 
@@ -207 +131 @@
     //}
 
-    Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
-    if ( Theta < ThetaMin )
+    Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
+    if ( theta < fThetaMin )
     {
       *fLog << "MSelBasic::Process; Run, Event, Theta = " 
             << fRawRun->GetRunNumber()<< ",  " 
-            << fMcEvt->GetEvtNumber() << ",  " << Theta << endl;
+            << fMcEvt->GetEvtNumber() << ",  " << theta << endl;
       rc = 1;
     }    
 
-    else if ( Theta > ThetaMax )
+    else if ( theta > fThetaMax )
     {
       rc = 2;
@@ -227 +151 @@
     }    
 
-    fErrors[rc]++;
+    fCut[rc]++;
 
     return rc==0 ? kTRUE : kCONTINUE;
@@ -246 +170 @@
     const Int_t entries = fEvt->GetNumPixels();
  
-    //$$$$$$$$$$$$$$$$$$
-    //const Int_t nall = fPed->GetSize();  
-    //*fLog << "nall = " << nall << endl;
-    //for (Int_t id=0; id<nall; id++)
-    //{
-    //  MGeomPix &gpix = (*fCam)[id];
-    //  if ( gpix.IsInOutermostRing() ) 
-    //  {
-    //    *fLog << "IsInOutermostRing : pixel no. = " << id << endl;
-    //  }
-
-    //  if ( gpix.IsInOuterRing() ) 
-    //  {
-    //    *fLog << "IsInOuterRing : pixel no. = " << id << endl;
-    //  }
-    //}
-    //$$$$$$$$$$$$$$$$$$
-
-
     for (Int_t i=0; i<entries; i++)
     {
@@ -316 +221 @@
     *fLog << GetDescriptor() << " execution statistics:" << endl;
     *fLog << dec << setfill(' ');
-    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
-          << (int)(fErrors[1]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: Zenith angle < " << ThetaMin << endl;
-
-    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) 
-          << (int)(fErrors[2]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: Zenith angle > " << ThetaMax << endl;
-
-    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) 
-          << (int)(fErrors[3]*100/GetNumExecutions()) 
+    *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) 
+          << (int)(fCut[1]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: Zenith angle < " << fThetaMin << endl;
+
+    *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) 
+          << (int)(fCut[2]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: Zenith angle > " << fThetaMax << endl;
+
+    *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) 
+          << (int)(fCut[3]*100/GetNumExecutions()) 
           << "%) Evts skipped due to: Software trigger not fullfilled" << endl;
 
-    *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) 
+    *fLog << " " << fCut[0] << " (" << (int)(fCut[0]*100/GetNumExecutions()) 
           << "%) Evts survived Basic selections!" << endl;
     *fLog << endl;
@@ -336 +241 @@
 
 
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MSelBasic.h

@@ -24 +24 @@
 {
 private:
-    const MPedestalCam *fPed;      // Pedestal information
-    const MGeomCam     *fCam;      // Camera Geometry 
-    const MCerPhotEvt  *fEvt;      // Cerenkov Photon Event 
-    const MMcEvt       *fMcEvt;       
+    const MPedestalCam  *fPed;      // Pedestal information
+    const MGeomCam      *fCam;      // Camera Geometry 
+    const MCerPhotEvt   *fEvt;      // Cerenkov Photon Event 
+    const MMcEvt        *fMcEvt;       
     const MRawRunHeader *fRawRun;       
 
-    Double_t     ThetaMin;
-    Double_t     ThetaMax;
+    Float_t     fThetaMin;
+    Float_t     fThetaMax;
 
-    Int_t        fErrors[4];
+    Int_t        fCut[4];
 
 public:
@@ -43 +43 @@
 
     Bool_t SwTrigger();
+    void SetCuts(Float_t thetamin, Float_t thetamax)
+         { fThetaMin = thetamin; fThetaMax = thetamax; }
 
     ClassDef(MSelBasic, 0)   // Task to evaluate basic cuts

trunk/MagicSoft/Mars/manalysis/MSelFinal.cc

@@ -67 +67 @@
     fHilSrcName = HilSrcName;
 
-    fHadronnessCut =   0.2;
-    fAlphaCut      = 100.0; //degrees
+    // default values of cuts
+    fHadronnessMax =   1.0;
+    fAlphaMax      = 100.0; //degrees
 }
 
@@ -109 +110 @@
     }
 
-    memset(fErrors, 0, sizeof(fErrors));
+    memset(fCut, 0, sizeof(fCut));
 
     return kTRUE;
@@ -133 +134 @@
     Double_t h = fHadronness->GetHadronness();
 
-    if ( h>fHadronnessCut )
+    if ( h>fHadronnessMax )
     {
       //*fLog << "MSelFinal::Process; h, alpha = " << h << ",  " 
@@ -140 +141 @@
     }    
 
-    else if ( modalpha > fAlphaCut )
+    else if ( modalpha > fAlphaMax )
     {
       //*fLog << "MSelFinal::Process; h, alpha = " << h << ",  " 
@@ -147 +148 @@
     }    
 
-    fErrors[rc]++;
+    fCut[rc]++;
 
     return rc==0 ? kTRUE : kCONTINUE;
@@ -164 +165 @@
     *fLog << GetDescriptor() << " execution statistics:" << endl;
     *fLog << dec << setfill(' ');
-    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
-          << (int)(fErrors[1]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: g/h separation cut (" << fHadronnessCut
-          << ")" << endl;
+    *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) 
+          << (int)(fCut[1]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: hadroness > "<< fHadronnessMax << endl;
 
-    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) 
-          << (int)(fErrors[2]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: cut in ALPHA (" << fAlphaCut
-          << " degrees)" << endl;
+    *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) 
+          << (int)(fCut[2]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: |ALPHA| > " << fAlphaMax
+          << " [degrees]" << endl;
 
-    *fLog << " " << fErrors[0] << " (" 
-          << (int)(fErrors[0]*100/GetNumExecutions()) 
+    *fLog << " " << fCut[0] << " (" 
+          << (int)(fCut[0]*100/GetNumExecutions()) 
           << "%) Evts survived Final selections!" << endl;
     *fLog << endl;

trunk/MagicSoft/Mars/manalysis/MSelFinal.h

@@ -32 +32 @@
 
     Double_t     fMm2Deg;   // conversion mm to degrees in camera
-    Int_t        fErrors[3];
+    Int_t        fCut[3];
     TString      fHilName;
     TString      fHilSrcName;
  
-    Float_t      fHadronnessCut;
-    Float_t      fAlphaCut;
+    Float_t      fHadronnessMax;
+    Float_t      fAlphaMax;
 
 public:
@@ -47 +47 @@
     Bool_t PostProcess();
 
-    void SetHadronnessCut(Float_t hadcut) { fHadronnessCut = hadcut; }
-    void SetAlphaCut(Float_t alpha)       { fAlphaCut      = alpha;  }
+    void SetCuts(Float_t hadmax, Float_t alphamax) 
+         { fHadronnessMax = hadmax; fAlphaMax = alphamax;  }
 
     ClassDef(MSelFinal, 0)   // Task to evaluate final cuts

trunk/MagicSoft/Mars/manalysis/MSelStandard.cc

@@ -64 +64 @@
     fHilName    = HilName;
     fHilSrcName = HilSrcName;
+
+    // default values of cuts
+    fUsedPixelsMax =   92;
+    fCorePixelsMin =    4;
+    fSizeMin       =   60;
+    fDistMin       =  0.4;
+    fDistMax       = 1.05;
+    fLengthMin     =  0.0;
+    fWidthMin      =  0.0;
 }
 
@@ -114 +123 @@
     //*fLog << "fMm2Deg = " << fMm2Deg << endl;
 
-    memset(fErrors, 0, sizeof(fErrors));
+    memset(fCut, 0, sizeof(fCut));
 
     return kTRUE;
@@ -131 +140 @@
     Int_t rc = 0;
 
-    Double_t fLength       = fHil->GetLength() * fMm2Deg;
-    Double_t fWidth        = fHil->GetWidth()  * fMm2Deg;
-    Double_t fDist         = fHilSrc->GetDist()* fMm2Deg;
-    //Double_t fDelta        = fHil->GetDelta()  * kRad2Deg;
-    Double_t fSize         = fHil->GetSize();
-    Int_t fNumUsedPixels   = fHil->GetNumUsedPixels();
-    Int_t fNumCorePixels   = fHil->GetNumCorePixels();
-
-    if ( fNumUsedPixels >= 92  ||  fNumCorePixels <= 4 )
-    {
-      //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = "
-      //      << fSize << ",  " << fDist << ",  " << fNumUsedPixels << ",  "
-      //      << fNumCorePixels << endl;
+    Double_t length       = fHil->GetLength() * fMm2Deg;
+    Double_t width        = fHil->GetWidth()  * fMm2Deg;
+    Double_t dist         = fHilSrc->GetDist()* fMm2Deg;
+    //Double_t delta        = fHil->GetDelta()  * kRad2Deg;
+    Double_t size         = fHil->GetSize();
+    Int_t numusedpixels   = fHil->GetNumUsedPixels();
+    Int_t numcorepixels   = fHil->GetNumCorePixels();
+
+    if ( numusedpixels >= fUsedPixelsMax  ||  numcorepixels <= fCorePixelsMin )
+    {
       rc = 1;
     }    
 
-    else if ( fSize <= 60.0         ||  fDist< 0.4           ||  fDist > 1.05 )
-    {
-      //*fLog << "MSelStandard::Process; fSize, fDist, fNumUsedPixels, fNumCorePixels = "
-      //      << fSize << ",  " << fDist << ",  " << fNumUsedPixels << ",  "
-      //      << fNumCorePixels << endl;
+    else if ( size <= fSizeMin )
+    {
       rc = 2;
     }    
 
-    else if ( fLength <= 0.0         ||  fWidth <= 0.0 )
-    {
-      //*fLog << "MSelStandard::Process; fLength, fWidth = "
-      //      << fLength << ",  " << fWidth << endl;
+    else if ( dist< fDistMin   ||  dist > fDistMax )
+    {
       rc = 3;
     }    
 
-
-    fErrors[rc]++;
+    else if ( length <= fLengthMin   ||  width <= fWidthMin )
+    {
+      rc = 4;
+    }    
+
+
+    fCut[rc]++;
 
     return rc==0 ? kTRUE : kCONTINUE;
@@ -180 +186 @@
     *fLog << GetDescriptor() << " execution statistics:" << endl;
     *fLog << dec << setfill(' ');
-    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
-          << (int)(fErrors[1]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: Requirements on no.of used or core pxels not fullfilled" << endl;
-
-    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) 
-          << (int)(fErrors[2]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: Requirements on SIZE or DIST not fullfilled" << endl;
-
-    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) 
-          << (int)(fErrors[3]*100/GetNumExecutions()) 
-          << "%) Evts skipped due to: Length or Width is <= 0" << endl;
-
-    *fLog << " " << fErrors[0] << " (" 
-          << (int)(fErrors[0]*100/GetNumExecutions()) 
+    *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) 
+          << (int)(fCut[1]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: Used pixels >= " << fUsedPixelsMax 
+          << " or Core pixels <= " << fCorePixelsMin << endl;
+
+    *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) 
+          << (int)(fCut[2]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: SIZE <= " << fSizeMin << endl;
+
+    *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) 
+          << (int)(fCut[3]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: DIST < " << fDistMin 
+          << " or DIST > " << fDistMax << endl;
+
+    *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) 
+          << (int)(fCut[4]*100/GetNumExecutions()) 
+          << "%) Evts skipped due to: LENGTH <= " << fLengthMin 
+          << " or WIDTH <= " << fWidthMin << endl;
+
+    *fLog << " " << fCut[0] << " (" 
+          << (int)(fCut[0]*100/GetNumExecutions()) 
           << "%) Evts survived Standard selections!" << endl;
     *fLog << endl;
@@ -201 +214 @@
 
 
+
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MSelStandard.h

@@ -30 +30 @@
 
     Double_t     fMm2Deg;   // conversion mm to degrees in camera
-    Int_t        fErrors[4];
+    Int_t        fCut[5];
     TString      fHilName;
     TString      fHilSrcName; 
+
+    Float_t     fUsedPixelsMax;
+    Float_t     fCorePixelsMin;
+    Float_t     fSizeMin;
+    Float_t     fDistMin;
+    Float_t     fDistMax;
+    Float_t     fLengthMin;
+    Float_t     fWidthMin;
 
 public:
@@ -42 +50 @@
     Bool_t PostProcess();
 
+    void SetCuts(Float_t usedpixelsmax, Float_t corepixelsmin,
+                 Float_t sizemin, Float_t distmin, Float_t distmax,
+                 Float_t lengthmin, Float_t widthmin)
+      { fUsedPixelsMax = usedpixelsmax; fCorePixelsMin = corepixelsmin;
+        fSizeMin = sizemin; fDistMin = distmin; fDistMax = distmax;
+        fLengthMin = lengthmin; fWidthMin = widthmin; } 
 
     ClassDef(MSelStandard, 0)   // Task to evaluate standard cuts

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 */
                  );
 

trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h

@@ -9 +9 @@
 #include "MRead.h"
 #endif
+
+#include <TRandom3.h>
 
 class TList;
@@ -22 +24 @@
 class MTaskList;
 class MParList;
+class MParameterD;
 
 struct outputpars;
@@ -42 +45 @@
     MRawRunHeader *fRawRunHeader; // raw run header
     MParList      *fParList;      // parameter list
+    MParameterD   *fHourAngle;    // hour angle [deg]
+    MParameterD   *fThetaOrig;    // original zenith angle [rad]
 
     Bool_t fIsMcFile;       // Flag whether current run is a MC run
@@ -53 +58 @@
     TArrayF fPedRMS;
 
+    TRandom3 ran3;
 
     Bool_t OpenNextFile();
@@ -77 +83 @@
     ~MCT1ReadPreProc();
 
-    Int_t AddFile(const char *fname, Int_t dummy=-1);
+    void AddFile(const char *fname);
 
     UInt_t GetEntries() { return fEntries; }
+
+    Bool_t SmearTheta(MParList *plist, Float_t *theta, Float_t *thetasmeared);
 
     ClassDef(MCT1ReadPreProc, 0) // Reads the CT1 preproc data file
@@ -86 +94 @@
 #endif
 
+
+

trunk/MagicSoft/Mars/mfileio/Makefile

@@ -20 +20 @@
 # @endcode 
 
-INCLUDES = -I. -I../mbase -I../mraw -I../mmc -I../mdata -I../manalysis -I../mgeom
+INCLUDES = -I. -I../mbase -I../mraw -I../mmc -I../mdata -I../manalysis -I../mgeom -I../mhist
 
 # @code