Changeset 1783 for trunk/MagicSoft/Mars

Timestamp:

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

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• macros/trigrate.C (modified) (2 diffs)
• mhist/MHMcRate.cc (modified) (4 diffs)
• mhist/MHMcRate.h (modified) (3 diffs)
• mmain/MMonteCarlo.cc (modified) (1 diff)
• mmontecarlo/MMcTriggerRateCalc.cc (modified) (8 diffs)
• mmontecarlo/MMcTriggerRateCalc.h (modified) (4 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
 
                                                          -*-*- END -*-*-
+ 2003/02/21: Abelardo Moralejo
+
+    * mmontecarlo/MMcTriggerRateCalc.[h,cc]
+      - adapted to new camera files, added warnings.
+      - added ReInit procedure to read relevant info from from the run headers
+
+    * mhist/MHMcRate.[h,cc]
+      - adapted accordingly. Added Set functions for several members.
+
+    * mmc/MMcCorsikaRunHeader.h
+      - added Get functions for fELowLim, fEUppLim and fSlopeSpec.
+
+    * mmain/MMontecarlo.cc, macros/trigrate.C
+      - adapted to changes above, changed MReadTree to MReadMarsFile to
+        be able to read the run headers.
+
  2003/02/21: Antonio Stamerra 
 

trunk/MagicSoft/Mars/macros/trigrate.C

@@ -124 +124 @@
     //    analised trigger conditions should be set (BgR[])
     //
-    MReadTree reader("Events", filename);
+
+    MReadMarsFile reader("Events", filename);
     tasklist.AddToList(&reader);
 
@@ -141 +142 @@
     cout << "Number of Trigger conditions: " << num << endl;
 
-    MMcTriggerRateCalc rate(dim, kPROTON, BgR, numnsbevents);
+    MMcTriggerRateCalc rate(dim, BgR, numnsbevents);
+
     tasklist.AddToList(&rate);
 

trunk/MagicSoft/Mars/mhist/MHMcRate.cc

@@ -18 +18 @@
 !   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@uni-sw.gwdg.de>
 !   Author(s): Harald Kornmayer 1/2001
+!   Author(s): Abelardo Moralejo 2/2003
+!
+!   Explanations on the rate calculation can be found in 
+!   chapter 7 of the following diploma thesis: 
+!   http://www.pd.infn.it/magic/tesi2.ps.gz (in Italian)
 !
 !   Copyright: MAGIC Software Development, 2000-2001
@@ -38 +43 @@
     fPartId=0;               // Type of particle
 
-    fEnergyMax=0.0;          // Maximum Energy in GeV
-    fEnergyMin=1000000.0;    // Minimum Energy in GeV
-
+    fEnergyMax=0.0;          // Maximum Energy (TeV)
+    fEnergyMin=1000000.0;    // Minimum Energy (TeV)
+
+    fSolidAngle = -1.;       // Solid angle within which incident directions
+                             // are distributed
     fThetaMax=0.0;           // Maximum theta angle of run
     fThetaMin=370.0;         // Minimum theta angle of run
@@ -119 +126 @@
 // --------------------------------------------------------------------------
 //
-//  Set the information about trigger due only to the background conditions
+//  Set the information about trigger due only to the Night Sky Background:
 //
 void MHMcRate::SetBackground (Float_t showers, Float_t triggers)
@@ -184 +191 @@
 
     if (fShowerRate <= 0)
-        fShowerRate = fFlux0/specidx*(epowmin-epowmax);
+      fShowerRate = fFlux0/specidx*(epowmax-epowmin);
+
+    if (fSolidAngle < 0.)
+      fSolidAngle = (fPhiMax-fPhiMin)*(cos(fThetaMin)-cos(fThetaMax));
 
     if (fPartId!=1)
-        fShowerRate *= (fPhiMax-fPhiMin)*(cos(fThetaMax)-cos(fThetaMin));
-
-    const Double_t impactdiff = fImpactMax-fImpactMin;
-
-    fShowerRate *= TMath::Pi()*(impactdiff/100.0*impactdiff/100.0);
+      fShowerRate *= fSolidAngle;
+
+    fShowerRate *= TMath::Pi()*(fImpactMax/100.0*fImpactMax/100.0 - 
+                                fImpactMin/100.0*fImpactMin/100.0);
 
     fShowerRateError = sqrt(fShowerRate);

trunk/MagicSoft/Mars/mhist/MHMcRate.h

@@ -12 +12 @@
     UShort_t fPartId;           // Type of particle
 
-    Float_t fEnergyMax;         // Maximum Energy in GeV
-    Float_t fEnergyMin;         // Minimum Energy in GeV
+    Float_t fEnergyMax;         // Maximum Energy in TeV
+    Float_t fEnergyMin;         // Minimum Energy in TeV
 
     Float_t fThetaMax;          // Maximum theta angle of run
@@ -20 +20 @@
     Float_t fPhiMin;            // Minimum phi angle of run
 
-    Float_t fImpactMax;         // Maximum impact parameter
-    Float_t fImpactMin;         // Minimum impact parameter
+    Float_t fSolidAngle;        // Solid angle within which incident directions
+                                // are distributed (sr)
+
+    Float_t fImpactMax;         // Maximum impact parameter (cm)
+    Float_t fImpactMin;         // Minimum impact parameter (cm)
 
     Float_t fBackTrig;          // Number of triggers from background
@@ -50 +53 @@
     void SetIncidentRate(Float_t showerrate);
 
+    void SetImpactMax(Float_t Impact) {fImpactMax=Impact;}
+    void SetImpactMin(Float_t Impact) {fImpactMin=Impact;}
+
+    void SetThetaMax(Float_t Theta) {fThetaMax=Theta;}
+    void SetThetaMin(Float_t Theta) {fThetaMin=Theta;}
+    void SetPhiMax(Float_t Phi) {fPhiMax=Phi;}
+    void SetPhiMin(Float_t Phi) {fPhiMin=Phi;}
+
+    void SetSolidAngle(Float_t Solid) {fSolidAngle=Solid;}
+    void SetEnergyMax(Float_t Energy) {fEnergyMax=Energy;}
+    void SetEnergyMin(Float_t Energy) {fEnergyMin=Energy;}
+
     void UpdateBoundaries(Float_t energy, Float_t theta, Float_t phi, Float_t impact);
 

trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc

@@ -349 +349 @@
     //    analised trigger conditions should be set (BgR[])
     //
-    MReadTree read("Events", fInputFile);
+    MReadMarsFile read("Events", fInputFile);
     tlist.AddToList(&read);
 
-    Float_t BgR[10]={660, 4, 0, 0, 0, 0, 0, 0, 0, 0};
-
-    MMcTriggerRateCalc crate(dim, 14, BgR, 100000);
+    // We calculate only shower rate (not including NSB-only triggers)
+    Float_t* BgR = new Float_t[dim];
+    memset(BgR, 0, num*sizeof(Float_t));
+
+    MMcTriggerRateCalc crate(dim, BgR, 100000);
     tlist.AddToList(&crate);
 

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

trunk/MagicSoft/Mars/mmontecarlo/MMcTriggerRateCalc.h

@@ -11 +11 @@
 class MParList;
 class MMcEvt;
+class MMcRunHeader;
+class MMcCorsikaRunHeader;
 class MMcTrig;
 class MHMcRate;
@@ -22 +24 @@
     TObjArray *fMcTrig;
 
+    MMcRunHeader *fMcRunHeader; 
+    MMcCorsikaRunHeader *fMcCorRunHeader; 
+
     UInt_t     fNum;           // decoded dimension
     UInt_t     fFirst;
     UInt_t     fLast;
 
-    Float_t*    fTrigger;   // Number of triggered showers
+    Float_t*   fTrigNSB;   // Number of triggers due to NSB alone
+    Float_t    fSimNSB;    // Number of simulated NSB-only events
 
+    Float_t*   fTrigger;       // Number of triggered showers
     Float_t    fShowers;       // Number of simulated showers
     Float_t    fAnalShow;      // Number of analysed showers
@@ -33 +40 @@
     Int_t      fPartId;        // Incident particle that generates showers
 
-    void Init(int dim, int part, float *trigbg,
+    void Init(int dim, float *trigbg,
               float simbg, const char *name, const char *title);
 
@@ -40 +47 @@
 
 public:
-    MMcTriggerRateCalc(int dim=0, int part=14, float *trigbg=NULL,
-                       float simbg=100000,
+    MMcTriggerRateCalc(int dim=0, float *trigbg=NULL, float simbg=100000,
                        const char *name=NULL, const char *title=NULL);
 
-    MMcTriggerRateCalc(float rate, int dim, int part, float *trigbg,
-                       float simbg,
+    MMcTriggerRateCalc(float rate, int dim, float *trigbg, float simbg,
                        const char *name=NULL, const char *title=NULL);
 
     ~MMcTriggerRateCalc();
+
+    Bool_t ReInit(MParList *plist);
 
     Bool_t PreProcess(MParList *pList);