Changeset 2036

Timestamp:

04/29/03 09:18:03 (21 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• mhistmc/MHMcCollectionArea.cc (modified) (6 diffs)
• mhistmc/MHMcCollectionArea.h (modified) (1 diff)
• mmontecarlo/MMcCollectionAreaCalc.cc (modified) (6 diffs)
• mmontecarlo/MMcCollectionAreaCalc.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+ 2003/04/29: Abelardo Moralej0
+
+   * mmontecarlo/MMcCollectionareaCalc.[h,cc]
+     - Now the energy limits and spectral slope used in Corsika are
+       read from the run header. Warning: impact parameter is still
+       fixed (I changed it from 400 m to 300 m, our present MC 
+       production default). Added warning message about this. 
+       Collection areas calculated up to now with the new MC library 
+       are overestimated by a ~77% ! This class will be removed in 
+       the future and replaced by MFillH.
+
+   * mhistmc/MHMcCollectionArea.[h,cc]
+     - Adapted to the change of MMcCollectionAreaCalc. Added arguments
+       to CalcEfficiency.
+
  2003/04/28: Nadia Tonello 
 
-  * mimage/MImgCleanStd.[h,cc]
+   * mimage/MImgCleanStd.[h,cc]
      - added the option kDemocratic using sigmabar of the 
         inner pixels
@@ -9 +24 @@
      - added documentation
 
-  * manalysis/MCerPhotPix.[h,cc]
+   * manalysis/MCerPhotPix.[h,cc]
      - added fRing and Get-Set functions 
  
  2003/04/28: Oscar Blanch
 
-  * mgeom/MGeomCamECO1000.[h,cc]
-    - New files for ECO 1000 camera geometry.
-
-  * mgeom/GeomLinkDef.h
-    - Added class MGeomCamECO1000
-
-  * mgeom/Makefile 
-    - Added class MGeomCamECO1000
+   * mgeom/MGeomCamECO1000.[h,cc]
+     - New files for ECO 1000 camera geometry.
+
+   * mgeom/GeomLinkDef.h
+     - Added class MGeomCamECO1000
+
+   * mgeom/Makefile 
+     - Added class MGeomCamECO1000
 
  2003/04/28: Abelardo Moralejo

trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.cc

@@ -40 +40 @@
 #include "MHMcEnergyImpact.h"
 
+#include "MLog.h"
+#include "MLogManip.h"
+
 ClassImp(MHMcCollectionArea);
 
@@ -223 +226 @@
 //  flag
 //
-void MHMcCollectionArea::CalcEfficiency(UInt_t numevts, Float_t angle)
+void MHMcCollectionArea::CalcEfficiency(UInt_t numevts, Float_t emin, Float_t emax, Float_t index)
 {
     // Here we estimate the total number of showers in each energy bin
@@ -238 +241 @@
     MH::SetBinning(fHistCol, &xaxis);
 
-    // TO BE FIXED! In forthcoming camera versions emin, emax and index
-    // will be available in a run header, and should be read from it.
-
-    const Float_t emin  = 10.;
-    const Float_t emax  = 30000.;  // Energies in GeV.
-    const Float_t index = 2.6;     // Differential spectral Index.
-
-    const Float_t expo = 1.-index;
+    const Float_t expo = 1+index;
 
     const Float_t k = (Float_t)numevts / (pow(emax,expo) - pow(emin,expo));
@@ -267 +263 @@
 
     // -----------------------------------------------------------
-
-    // Impact parameter range:
+    //
+    // Impact parameter range:  TO BE FIXED! Impact parameter shoud be
+    // read from run header, but it is not yet in!!
+    //
     const Float_t r1 = 0;
-    const Float_t r2 = 400;
-
-    const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);
-
-    angle *= TMath::Pi()/180;
+    const Float_t r2 = 300;
+
+    *fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of 300 meters for the MC events. Check that this is the true one!" <<endl<<endl;
+    const Float_t area = TMath::Pi() * (r2*r2 - r1*r1);
 
     for (Int_t ix=1; ix<=nbinx; ix++)
@@ -294 +291 @@
         const Double_t eff = Ns/Na;
 
-        const Double_t err = sqrt((1.-eff)*Ns)/Na;
-
-        const Float_t area = dr * cos(angle);
+        const Double_t efferr = sqrt((1.-eff)*Ns)/Na;
 
         fHistCol->SetBinContent(ix, eff*area);
-        fHistCol->SetBinError(ix, err*area);
+        fHistCol->SetBinError(ix, efferr*area);
     }
 
@@ -340 +335 @@
 
             const Double_t eff = h.GetCellContent(ix, iy);
-            const Double_t err = h.GetCellError(ix, iy);
+            const Double_t efferr = h.GetCellError(ix, iy);
 
             colA += eff*A;
-            errA += A*A * err*err;
+            errA += A*A * efferr*efferr;
         }
 

trunk/MagicSoft/Mars/mhistmc/MHMcCollectionArea.h

@@ -38 +38 @@
 
     void CalcEfficiency();
-    void CalcEfficiency(UInt_t allevts, Float_t theta);
+    void CalcEfficiency(UInt_t allevts, Float_t emin, Float_t emax, Float_t index);
 
     void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall);

trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc

@@ -28 +28 @@
 //  MHMcCollectionAreaCalc
 //
-//  Remark: The initialization is maily done in the ReInit function.
+//  Remark: The initialization is mainly done in the ReInit function.
 //          Please make sure, that you don't use MReadTree when processing
 //          a file. Use a 'ReInit'-calling task like MReadMarsFile
@@ -44 +44 @@
 #include "MMcTrig.hxx"
 #include "MMcRunHeader.hxx"
+#include "MMcCorsikaRunHeader.h"
 
 #include "MHMcCollectionArea.h"
@@ -77 +78 @@
         return kFALSE;
 
-    fTheta                    = -1;
+    fTheta                    =  -1;
+    fEmin                     =  -1;
+    fEmax                     =  -1;
+    fSlope                    =   0;
     fTotalNumSimulatedShowers =  0;
     fCorsikaVersion           =  0;
@@ -91 +95 @@
     {
         *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl;
+        return kFALSE;
+    }
+
+    MMcCorsikaRunHeader *corrunheader  = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
+    if (!corrunheader)
+    {
+        *fLog << err << dbginf << "Error - MMcCorsikaRunHeader not found... exit." << endl;
         return kFALSE;
     }
@@ -110 +121 @@
 
     fCorsikaVersion = runheader->GetCorsikaVersion();
+
+    if ( fEmin > 0 &&
+         (fEmin  != corrunheader->GetELowLim() ||
+          fEmax  != corrunheader->GetEUppLim() ||
+          fSlope != corrunheader->GetSlopeSpec()) )
+      *fLog << warn << dbginf << "Warning - Read files have different energy distribution..." << endl;
+
+    fEmin  = corrunheader->GetELowLim();
+    fEmax  = corrunheader->GetEUppLim();
+    fSlope = corrunheader->GetSlopeSpec();
 
     fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
@@ -167 +188 @@
     {
         *fLog << inf << "Total number of showers: " << fTotalNumSimulatedShowers << endl;
-        fCollArea->CalcEfficiency(fTotalNumSimulatedShowers,
-                                  fCorsikaVersion == 5200 ? fTheta : 0);
+        fCollArea->CalcEfficiency(fTotalNumSimulatedShowers, fEmin, fEmax, fSlope);
     }
 

trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h

@@ -24 +24 @@
     UInt_t fCorsikaVersion;
     Float_t fTheta;
+    Float_t fEmin;
+    Float_t fEmax;
+    Float_t fSlope;
 
     Bool_t fAllEvtsTriggered;