Changeset 19535 for trunk/Mars/mcorsika

Timestamp:

06/03/19 15:55:18 (5 years ago)

Author:

tbretz

Message:

Allow to read telescope positions from a corsika input card and apply them to plain corsika data according to the telescope is in the filename.

Location:

trunk/Mars/mcorsika

Files:

• MCorsikaRead.cc (modified) (7 diffs)
• MCorsikaRead.h (modified) (1 diff)

trunk/Mars/mcorsika/MCorsikaRead.cc

@@ -43 +43 @@
 
 #include <TSystem.h>
+#include <TPRegexp.h>
+#include <TObjString.h>
 
 #include "MLog.h"
@@ -100 +102 @@
     : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fArrayIdx(-1), 
     fForceMode(kFALSE), fFileNames(0), fNumFile(0), fNumEvents(0), 
-    fNumTotalEvents(0), fInFormat(0), fParList(0), fNumTelescopes(1), fReadState(0)
+    fNumTotalEvents(0), fInFormat(0), fParList(0), fNumTelescopes(1), fNumTelescope(0), fReadState(0)
 {
     fName  = name  ? name  : "MRead";
@@ -195 +197 @@
 
     *fLog << inf << "Open file: '" << name << "'" << endl;
+
+    // Check if thsi is a single telescope from a CORSIKA telescope array
+    // This is indicated by the -telXXX attached to the filename
+    TObjArray *res = TPRegexp("(.*/)?cer([0-9]+)-tel([0-9]+)").MatchS(name, "i");
+    if (res->GetLast()==3)
+    {
+        const TString match = res->At(0)->GetName();
+        const TString path  = res->At(1)->GetName();
+        const UInt_t  runid = atoi(res->At(2)->GetName());
+        const UInt_t  telid = atoi(res->At(3)->GetName());
+
+        *fLog << inf2 << "Run " << runid << " detected to be for the " << telid << "-th telescope." << endl;
+
+        *fLog << inf << "Telescope #" << telid << " [X/Y/Z (R)]: ";
+        *fLog << fTelescopeX[telid-1] << "/";
+        *fLog << fTelescopeY[telid-1] << "/";
+        *fLog << fTelescopeZ[telid-1] << " (R=" << fTelescopeR[telid-1] << ")" << endl;
+
+        fNumTelescope = telid;
+    }
+    delete res;
+
 
     if (fDisplay)
@@ -583 +607 @@
                   Float_t xArrOff, yArrOff;
                   fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
-                  fEvtHeader->SetTelescopeOffset(arrayIdx, 
+                  fEvtHeader->SetTelescopeOffset(arrayIdx,
                                                    xArrOff + fTelescopeY[telIdx], 
                                                    yArrOff - fTelescopeX[telIdx] );
@@ -605 +629 @@
                   fEvtHeader->ResetNumReuse();
                   fEvtHeader->InitXY();
+                  if (fNumTelescope>0)
+                  {
+                      // Here, the impact has opposite sign -- I don't understand why
+                      fEvtHeader->SetXY(+fTelescopeY[fNumTelescope-1]-fEvtHeader->GetX(),
+                                        -fTelescopeX[fNumTelescope-1]-fEvtHeader->GetY());
+                  }
                   fBlockType = 1109;  // save corsika events
                   continue;
@@ -640 +670 @@
 
             case 1109:  // save corsika events
-               fEvtHeader->InitXY();
+                fEvtHeader->InitXY();
+                if (fNumTelescope>0)
+                {
+                      // Here, the impact has opposite sign -- I don't understand why
+                      fEvtHeader->SetXY(+fTelescopeY[fNumTelescope-1]-fEvtHeader->GetX(),
+                                        -fTelescopeX[fNumTelescope-1]-fEvtHeader->GetY());
+                }
                status = fEvent->ReadCorsikaEvt(fRawEvemtBuffer,
                                                fBlockLength == MCorsikaFormat::kBlockLengthRaw/21 - 4 ? 7 : 8,
@@ -726 +762 @@
     //*fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
 
-
     return kTRUE;
 }
 
+Int_t MCorsikaRead::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "TelescopeIndex", print))
+    {
+        rc = kTRUE;
+        fNumTelescope = GetEnvValue(env, prefix, "TelescopeIndex", 0);
+
+        fTelescopeX.Set(1);
+        fTelescopeY.Set(1);
+        fNumTelescope = 1;
+
+        if (IsEnvDefined(env, prefix, "TelescopeX", print))
+            fTelescopeX[0] = GetEnvValue(env, prefix, "TelescopeX", 0.);
+        if (IsEnvDefined(env, prefix, "TelescopeY", print))
+            fTelescopeY[0] = GetEnvValue(env, prefix, "TelescopeY", 0.);
+
+        *fLog << all << "Telescope #" << setw(4) << fNumTelescope << " [X/Y]: ";
+        *fLog << setw(7) << fTelescopeX[0] << "/";
+        *fLog << setw(7) << fTelescopeY[0] << endl;
+    }
+
+    // Read telescope positions from corsika input card
+    if (IsEnvDefined(env, prefix, "CorsikaInputCard", print))
+    {
+        rc = kTRUE;
+        TString fname = GetEnvValue(env, prefix, "CorsikaInputCard", "");
+
+        gSystem->ExpandPathName(fname);
+
+        // Is file accessible
+        if (gSystem->AccessPathName(fname, kFileExists))
+            return kERROR;
+
+        *fLog << inf2 << "Reading Telecope positions from " << fname << ":" << endl;
+
+        TPRegexp regexp("\\s*TELESCOPE"
+                        "\\s+([-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?)"
+                        "\\s+([-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?)"
+                        "\\s+([-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?)"
+                        "\\s+([-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?)"
+                        "\\s+([0-9]+)\\s*");
+
+        ifstream fin(fname);
+
+
+        TString buf;
+        while (1)
+        {
+            buf.ReadLine(fin);
+            if (!fin)
+                break;
+
+            const TObjArray *res = regexp.MatchS(buf);
+            if (res->GetLast()==9)
+            {
+                const UInt_t idx = atoi(res->At(9)->GetName());
+                if (idx==0)
+                {
+                    *fLog << warn << "Telescope with index 0 ignored." << endl;
+                    delete res;
+                    continue;
+                }
+                if (idx>=fNumTelescopes)
+                {
+                    fTelescopeX.Set(idx);
+                    fTelescopeY.Set(idx);
+                    fTelescopeZ.Set(idx);
+                    fTelescopeR.Set(idx);
+                    fNumTelescopes = idx;
+                }
+
+                fTelescopeX[idx-1] = atof(res->At(1)->GetName());
+                fTelescopeY[idx-1] = atof(res->At(3)->GetName());
+                fTelescopeZ[idx-1] = atof(res->At(5)->GetName());
+                fTelescopeR[idx-1] = atof(res->At(7)->GetName());
+            }
+            delete res;
+        }
+
+
+        for (int i=0; i<fNumTelescopes; i++)
+        {
+            *fLog << all << "Telescope #" << setw(4) << i+1 << " [X/Y/Z (R)]: ";
+            *fLog << setw(7) << fTelescopeX[i] << "/";
+            *fLog << setw(7) << fTelescopeY[i] << "/";
+            *fLog << setw(7) << fTelescopeZ[i] << "  (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
+        }
+    }
+
+    return rc;
+}
+

trunk/Mars/mcorsika/MCorsikaRead.h

@@ -93 +93 @@
     UInt_t GetEntries() { return fNumTotalEvents/*/fInterleave*/; }
 
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
+
     ClassDef(MCorsikaRead, 0)   // Task to read the raw data binary file
 };