Changeset 1783 for trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc

Timestamp:

02/21/03 18:09:36 (22 years ago)

Author:

moralejo

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc (modified) (8 diffs)

trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.cc

@@ -38 +38 @@
 #include "MMcEvt.hxx"
 #include "MMcTrig.hxx"
+#include "MMcRunHeader.hxx"
+#include "MMcCorsikaRunHeader.h"
 
 ClassImp(MMcTriggerRateCalc);
 
-void MMcTriggerRateCalc::Init(int dim, int part, float *trigbg,
-                              float simbg,
+void MMcTriggerRateCalc::Init(int dim, float *trigbg, float simbg,
                               const char *name, const char *title)
 {
@@ -51 +52 @@
     fMcRate = NULL;
 
+    fTrigNSB = trigbg;
+    fSimNSB = simbg;
+
+    fPartId = -1; 
+
     fShowers = 0;
-    fAnalShow = simbg;
-
-    fPartId=part;
+    fAnalShow = 0;
 
     fFirst = dim>0 ?   1 : -dim;
@@ -60 +64 @@
 
     fNum = fLast-fFirst+1;
-
     fTrigger = new float[fNum];
 
     for (UInt_t i=0;i<fNum;i++)
-        fTrigger[i] = dim&&trigbg ? trigbg[i] : 0;
-
+      fTrigger[i]=0;
+
+    AddToBranchList("MMcEvt.fPartId");
     AddToBranchList("MMcEvt.fImpact");
     AddToBranchList("MMcEvt.fEnergy");
@@ -72 +76 @@
     AddToBranchList("MMcEvt.fPhotElfromShower");
     AddToBranchList("MMcTrig", "fNumFirstLevel", fFirst, fLast);
+
+}
+
+// ReInit: read run header to get some info later:
+
+Bool_t MMcTriggerRateCalc::ReInit(MParList *pList)
+{
+    fMcRunHeader = (MMcRunHeader*) pList->FindObject("MMcRunHeader");
+    if (!fMcRunHeader)
+    {
+        *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl;
+        return kFALSE;
+    }
+
+    fMcCorRunHeader = (MMcCorsikaRunHeader*) pList->FindObject("MMcCorsikaRunHeader");
+    if (!fMcCorRunHeader)
+    {
+        *fLog << err << dbginf << "Error - MMcCorsikaRunHeader not found... exit." << endl;
+        return kFALSE;
+    }
+
+    if (fMcRunHeader->GetAllEvtsTriggered())
+        fShowers += fMcRunHeader->GetNumSimulatedShowers();
+
+    if (fMcRunHeader->GetAllEvtsTriggered())
+      {
+        *fLog << endl << all << endl <<
+          "WARNING: the input data file contains only the" << endl << 
+          "events that triggered. I will assume the standard value" << endl <<
+          "for maximum impact parameter (400 m)" <<endl;
+
+
+        if (fTrigNSB[0] > 0)
+          *fLog << endl << all <<
+            "WARNING: NSB rate can be overestimated by up to 5%." << endl <<
+            "For a precise estimate of the total rate including NSB" << endl << 
+            "accidental triggers I need a file containing all event headers." 
+                << endl;
+        else 
+          *fLog << endl << all <<
+            "WARNING: calculating only shower rate (no NSB accidental triggers)" << endl;
+      }
+
+    *fLog << endl << all <<
+      "WARNING: I will assume the standard maximum off axis angle" << endl <<
+      "(5 degrees) for the original showers" << endl;
+
+    for (UInt_t i=0; i<fNum; i++)
+      {
+        if (fMcRunHeader->GetAllEvtsTriggered())
+          {
+            GetRate(i)->SetImpactMin(0.);
+            GetRate(i)->SetImpactMax(40000.);   // in cm
+          }
+        GetRate(i)->SetSolidAngle(2.390941702e-2);  // sr
+
+        //
+        // Set limits of energy range, coverting from GeV to TeV:
+        //
+        GetRate(i)->SetEnergyMin(fMcCorRunHeader->GetELowLim()/1000.);
+        GetRate(i)->SetEnergyMax(fMcCorRunHeader->GetEUppLim()/1000.);
+
+      }
+
+    return kTRUE;
 }
 
@@ -79 +148 @@
 //
 //      dim:     fDimension
-//      part:    fPartId
-//      *trigbg: number of shower from bacground that triggers
+//      *trigbg: number of NSB triggers for
 //               a given trigger condition.
-//      simbg:   Number of simulated showers for the background
+//      simbg:   Number of simulated events for the NSB
 //      rate:    rate of incident showers
 //
 
-MMcTriggerRateCalc::MMcTriggerRateCalc(float rate, int dim, int part,
+MMcTriggerRateCalc::MMcTriggerRateCalc(float rate, int dim,
                                        float *trigbg, float simbg,
                                        const char *name, const char *title)
 {
-    Init(dim, part, trigbg, simbg, name, title);
+    Init(dim, trigbg, simbg, name, title);
 }
 
@@ -99 +167 @@
 //
 //      dim:     fDimension
-//      part:    fPartId
-//      *trigbg: number of shower from bacground that triggers
+//      *trigbg: number of NSB triggers for
 //               a given trigger condition.
-//      simbg:   Number of simulated showers for the background
-//
-MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, int part, float *trigbg,
-                                       float simbg,
+//      simbg:   Number of simulated events for the NSB
+//
+MMcTriggerRateCalc::MMcTriggerRateCalc(int dim, float *trigbg,float simbg,
                                        const char *name, const char *title)
 {
-    Init(dim, part, trigbg, simbg, name, title);
+    Init(dim, trigbg, simbg, name, title);
 }
 
@@ -160 +226 @@
     }
 
+    for (UInt_t i=0;i<fNum;i++)
+      {
+        MHMcRate &rate = *GetRate(i);
+
+        if (fTrigNSB)
+          rate.SetBackground(fTrigNSB[i], fSimNSB);
+        else
+          rate.SetBackground(0., fSimNSB);
+      }
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+//  The Process-function counts the number of simulated showers, the
+//  number of analised showers and the number of triggers. It also updates
+//  the limits for theta, phi, energy and impact parameter in the
+//  MHMcRate container.
+//
+Bool_t MMcTriggerRateCalc::Process()
+{
+    //
+    //  Counting analysed showers (except in the case in which the file 
+    //  contains only events that triggered, since then we do not know
+    //  how many showers were analysed).
+    //
+
+    if ( ! fMcRunHeader->GetAllEvtsTriggered() &&
+         fMcEvt->GetPhotElfromShower() )
+      fAnalShow++;
+
+    //
+    // Set PartId and check it is the same for all events
+    //
+    if (fPartId < 0)
+      fPartId = fMcEvt->GetPartId();
+    else if (fPartId != fMcEvt->GetPartId())
+      {
+        *fLog << err << dbginf << "Incident particle type seems to change";
+        *fLog << "from " << fPartId << " to " << fMcEvt->GetPartId() << endl;
+        *fLog << "in input data files... aborting." << endl;
+        return kFALSE;
+      }
+
+    //
+    //  Getting angles, energy and impact parameter to set boundaries
+    //
+    const Float_t theta = fMcEvt->GetTheta();
+    const Float_t phi   = fMcEvt->GetPhi();
+    const Float_t param = fMcEvt->GetImpact();
+    const Float_t ener  = fMcEvt->GetEnergy()/1000.0;
+
+    //
+    //  Counting number of triggers
+    //
+    for (UInt_t i=0; i<fNum; i++)
+    {
+        if (GetTrig(i)->GetFirstLevel())
+            fTrigger[i] ++; 
+
+        if ( ! fMcRunHeader->GetAllEvtsTriggered())
+          GetRate(i)->UpdateBoundaries(ener, theta, phi, param);
+    }
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+//  The PostProcess-function calculates and shows the trigger rate
+//
+Bool_t MMcTriggerRateCalc::PostProcess()
+{
     for (UInt_t i=0; i<fNum; i++)
     {
@@ -168 +308 @@
         {
         case kPROTON:
-            rate.SetFlux(0.1091, 2.75);
-            break;
+          if (fMcCorRunHeader->GetSlopeSpec() != -2.75)
+            {
+              *fLog << err << dbginf <<
+                "Spectrum slope as read from  input file ("<< 
+                fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) for protons" << endl << "... aborting." << endl;
+              return kFALSE;
+            }
+          rate.SetFlux(0.1091, 2.75);
+          break;
         case kHELIUM:
-            rate.SetFlux(0.0660, 2.62);
-            break;
+          if (fMcCorRunHeader->GetSlopeSpec() != -2.62)
+            {
+              *fLog << err << dbginf <<
+                "Spectrum slope as read from  input file ("<< 
+                fMcCorRunHeader->GetSlopeSpec() << ") does not match the expected one (-2.75) for Helium" << endl << "... aborting." << endl;
+              return kFALSE;
+            }
+          rate.SetFlux(0.0660, 2.62);
+          break;
         default:
-            *fLog << err << dbginf << "Unknown incident flux parameters for ";
-            *fLog << fPartId<< " particle Id ... aborting." << endl;
-            return kFALSE;
+          *fLog << err << dbginf << "Unknown incident flux parameters for ";
+          *fLog << fPartId<< " particle Id ... aborting." << endl;
+          return kFALSE;
         }
-        rate.SetBackground(fTrigger[i], fAnalShow);
-
-        fTrigger[i]=0;
-    }
-
-    fAnalShow=0.0;
-
-    return kTRUE;
-}
-
-// --------------------------------------------------------------------------
-//
-//  The Process-function counts the number of simulated showers, the
-//  number of analised showers and the number of triggers. It also updates
-//  the limits for theta, phi, energy and impact parameter in the
-//  MHMcRate container.
-//
-Bool_t MMcTriggerRateCalc::Process()
-{
-    //
-    //  Counting analysed and simulated showers
-    //
-    fShowers++;
-    if (fMcEvt->GetPhotElfromShower())
-        fAnalShow++;
-
-    //
-    //  Getting angles, energy and impact parameter to set boundaries
-    //
-    const Float_t theta = fMcEvt->GetTheta();
-    const Float_t phi   = fMcEvt->GetPhi();
-    const Float_t param = fMcEvt->GetImpact();
-    const Float_t ener  = fMcEvt->GetEnergy()/1000.0;
-
-    //
-    //  Counting number of triggers
-    //
-    for (UInt_t i=0; i<fNum; i++)
-    {
-        if (GetTrig(i)->GetFirstLevel())
-            fTrigger[i] ++; 
-
-        GetRate(i)->UpdateBoundaries(ener, theta, phi, param);
-    }
-
-    return kTRUE;
-}
-
-// --------------------------------------------------------------------------
-//
-//  The PostProcess-function calculates and shows the trigger rate
-//
-Bool_t MMcTriggerRateCalc::PostProcess()
-{
+    }
+
     //
     // Computing trigger rate