Changeset 19848 for trunk/Mars/mcorsika/MCorsikaRead.cc

Timestamp:

11/01/19 14:18:10 (5 years ago)

Author:

tbretz

Message:

Make a public access to MCorsikaRead::CalcNumEvents possible. This allows to use it to check whether all files are valid. Added the new option '--check' to ceres for that. Copy '-tel' to the file names. This is more conevenient in reading them.

File:

• trunk/Mars/mcorsika/MCorsikaRead.cc (modified) (9 diffs)

trunk/Mars/mcorsika/MCorsikaRead.cc

@@ -167 +167 @@
 // This opens the next file in the list and deletes its name from the list.
 //
-Int_t MCorsikaRead::OpenNextFile(Bool_t print)
+Int_t MCorsikaRead::OpenNextFile(Bool_t print, Bool_t telrequired)
 {
 
@@ -196 +196 @@
         return kERROR;
 
-    *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)
+    *fLog << (telrequired?inf:all) << "Open file: '" << name << "'" << endl;
+
+    if (telrequired)
     {
-        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 telescope #" << telid << endl;
-
-        if (telid>=fNumTelescopes)
+        // 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)
         {
-            *fLog << err << "Position of telescope " << telid << " not defined." << endl;
-            return kERROR;
+            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 telescope #" << telid << endl;
+
+            if (telid>=fNumTelescopes)
+            {
+                *fLog << err << "Position of telescope " << telid << " not defined." << endl;
+                return kERROR;
+            }
+
+            *fLog << inf << "Telescope #" << telid << " [X/Y/Z (R)]: ";
+            *fLog << fTelescopeX[telid] << "/";
+            *fLog << fTelescopeY[telid] << "/";
+            *fLog << fTelescopeZ[telid] << " (R=" << fTelescopeR[telid] << ")" << endl;
+
+            fNumTelescope = telid+1;
         }
-
-        *fLog << inf << "Telescope #" << telid << " [X/Y/Z (R)]: ";
-        *fLog << fTelescopeX[telid] << "/";
-        *fLog << fTelescopeY[telid] << "/";
-        *fLog << fTelescopeZ[telid] << " (R=" << fTelescopeR[telid] << ")" << endl;
-
-        fNumTelescope = telid+1;
+        delete res;
     }
-    delete res;
-
 
     if (fDisplay)
@@ -288 +290 @@
     fNumTotalEvents = 0;
 
+    MCorsikaRunHeader runheader;
+
     Bool_t rc = kTRUE;
 
     while (1)
     {
-        switch (OpenNextFile(kFALSE))
+        switch (OpenNextFile(kFALSE, kFALSE))
         {
         case kFALSE:
@@ -310 +314 @@
                {
                if (blockType == 1200)
-                  status = fRunHeader->ReadEvt(fInFormat, blockLength / sizeof(Float_t));
+                  status = runheader.ReadEvt(fInFormat, blockLength / sizeof(Float_t));
 
                else if(blockType == 1201)
@@ -319 +323 @@
                   vector<Float_t> buffer(blockLength / sizeof(Float_t));
                   status = fInFormat->Read(buffer.data(), blockLength);
-                  status = fRunHeader->ReadEventHeader(buffer.data());
+                  status = runheader.ReadEventHeader(buffer.data());
                   break;
                   }
@@ -325 +329 @@
                   fInFormat->Seek(blockLength);
                }
-                  
+
             if (status != kTRUE)
                return status;
@@ -331 +335 @@
             if (!fInFormat->SeekEvtEnd())
             {
-               *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl;
+               *fLog << (fForceMode?warn:err) << "ERROR - RUNE section not found in file." << endl;
                if (!fForceMode)
                   return fForceMode ? kTRUE : kFALSE;
             }
 
-            if (!fRunHeader->ReadEvtEnd(fInFormat, kTRUE))
+            if (!runheader.ReadEvtEnd(fInFormat, kTRUE))
             {
-               *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl;
+               *fLog << (fForceMode?warn:err) << "ERROR - Reading RUNE section failed." << endl;
                if (!fForceMode)
                   return kFALSE;
             }
 
-            fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()*
+            fNumTotalEvents += runheader.GetNumEvents()*runheader.GetNumReuse()*
                 (fTelescopeIdx<0 && fNumTelescope==0 ? fNumTelescopes : 1);
             continue;
@@ -806 +810 @@
         *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);
 
@@ -823 +820 @@
                 break;
 
-            const TObjArray *res = regexp.MatchS(buf);
-            if (res->GetLast()==9)
+            buf = buf.Strip(TString::kBoth);
+            if (buf(0, 9)!="TELESCOPE")
+                continue;
+
+            float x, y, z, r;
+            int idx;
+            if (sscanf(buf.Data()+9, "%f %f %f %f %d", &x, &y, &z, &r, &idx)!=5)
+                continue;
+
+            if (idx>=fNumTelescopes)
             {
-                const UInt_t idx = atoi(res->At(9)->GetName());
-                if (idx>=fNumTelescopes)
-                {
-                    fNumTelescopes = idx+1;
-                    fTelescopeX.Set(fNumTelescopes);
-                    fTelescopeY.Set(fNumTelescopes);
-                    fTelescopeZ.Set(fNumTelescopes);
-                    fTelescopeR.Set(fNumTelescopes);
-                }
-
-                fTelescopeX[idx] = atof(res->At(1)->GetName());
-                fTelescopeY[idx] = atof(res->At(3)->GetName());
-                fTelescopeZ[idx] = atof(res->At(5)->GetName());
-                fTelescopeR[idx] = atof(res->At(7)->GetName());
+                fNumTelescopes = idx+1;
+                fTelescopeX.Set(fNumTelescopes);
+                fTelescopeY.Set(fNumTelescopes);
+                fTelescopeZ.Set(fNumTelescopes);
+                fTelescopeR.Set(fNumTelescopes);
             }
-            delete res;
+
+            fTelescopeX[idx] = x;
+            fTelescopeY[idx] = y;
+            fTelescopeZ[idx] = z;
+            fTelescopeR[idx] = r;
         }