source: trunk/Mars/mcorsika/MCorsikaRead.cc@ 20076

Last change on this file since 20076 was 20043, checked in by tbretz, 4 years ago
Copy the telescope number to the coriska run-header
File size: 26.8 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 11/2008 <mailto:tbretz@astro.uni-wuerzburg.de>
19 Reiner Rohlfs 11/2010 <mailto:Reiner.Rohlfs@unige.ch>
20!
21! Copyright: Software Development, 2000-2008
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MCorsikaRead
29//
30// Input Containers:
31// -/-
32//
33// Output Containers:
34// MCorsikaRunHeader
35// MCorsikaEvtHeader
36// MPhotonEvent
37//
38//////////////////////////////////////////////////////////////////////////////
39#include "MCorsikaRead.h"
40
41#include <errno.h>
42#include <fstream>
43
44#include <TSystem.h>
45#include <TPRegexp.h>
46#include <TObjString.h>
47
48#include "MLog.h"
49#include "MLogManip.h"
50
51#include "MParList.h"
52#include "MStatusDisplay.h"
53
54#include "MCorsikaFormat.h"
55#include "MCorsikaRunHeader.h"
56#include "MCorsikaEvtHeader.h"
57
58#include "MPhotonEvent.h"
59
60ClassImp(MCorsikaRead);
61
62using namespace std;
63
64
65/* ----------- please don't delete and don't care about (Thomas) ------------
66#define kBUFSZ 64 //1024*1024*64
67#include <iomanip.h>
68class bifstream : public istream, public streambuf
69{
70private:
71 char fBuffer[kBUFSZ]; //!
72 FILE *fd;
73
74 int sync()
75 {
76 memset(fBuffer, 0, kBUFSZ);
77 return 0;
78 }
79 int underflow()
80 {
81 int sz=fread(fBuffer, kBUFSZ, 1, fd);
82 //int sz=fread(fBuffer, 1, kBUFSZ, fd);
83 setg(fBuffer, fBuffer, fBuffer+kBUFSZ);
84
85 return sz==1 ? *(unsigned char*)fBuffer : EOF;//EOF;
86 //return sz==kBUFSZ ? *(unsigned char*)fBuffer : EOF;//EOF;
87 }
88public:
89 bifstream(const char *name) : istream(this)
90 {
91 fd = fopen(name, "rb");
92 setbuf(fBuffer, kBUFSZ);
93 }
94};
95*/
96
97// --------------------------------------------------------------------------
98//
99// Default constructor. It tries to open the given file.
100//
101MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title)
102 : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fArrayIdx(-1),
103 fForceMode(kFALSE), fFileNames(0), fNumFile(0), fNumEvents(0),
104 fNumTotalEvents(0), fInFormat(0), fParList(0), fNumTelescopes(1), fNumTelescope(0), fReadState(0)
105{
106 fName = name ? name : "MRead";
107 fTitle = title ? title : "Read task to read DAQ binary files";
108
109 fFileNames = new TList;
110 fFileNames->SetOwner();
111
112 if (fname!=NULL)
113 AddFile(fname);
114}
115
116// --------------------------------------------------------------------------
117//
118// Destructor. Delete input stream.
119//
120MCorsikaRead::~MCorsikaRead()
121{
122 delete fFileNames;
123 if (fInFormat)
124 delete fInFormat;
125}
126
127/*
128Byte_t MCorsikaRead::IsFileValid(const char *name)
129{
130 MZlib fin(name);
131 if (!fin)
132 return 0;
133
134 Byte_t c[4];
135 fin.read((char*)c, 4);
136 if (!fin)
137 return 0;
138
139 if (c[0]!=0xc0)
140 return 0;
141
142 if (c[1]==0xc0)
143 return 1;
144
145 if (c[1]==0xc1)
146 return 2;
147
148 return 0;
149}
150*/
151
152// --------------------------------------------------------------------------
153//
154// Add a new file to a list of files to be processed, Returns the number
155// of files added. (We can enhance this with a existance chack and
156// wildcard support)
157//
158Int_t MCorsikaRead::AddFile(const char *fname, Int_t entries)
159{
160 TNamed *name = new TNamed(fname, "");
161 fFileNames->AddLast(name);
162 return 1;
163
164}
165// --------------------------------------------------------------------------
166//
167// This opens the next file in the list and deletes its name from the list.
168//
169Int_t MCorsikaRead::OpenNextFile(Bool_t print, Bool_t telrequired)
170{
171
172 //
173 // open the input stream and check if it is really open (file exists?)
174 //
175 if (fInFormat)
176 delete fInFormat;
177 fInFormat = NULL;
178
179 //
180 // Check for the existance of a next file to read
181 //
182 TObject *file = fFileNames->At(fNumFile);
183 if (!file)
184 return kFALSE;
185
186 //
187 // open the file which is the first one in the chain
188 //
189 const char *name = file->GetName();
190
191 const char *expname = gSystem->ExpandPathName(name);
192 fInFormat = MCorsikaFormat::CorsikaFormatFactory(expname);
193 delete [] expname;
194
195 if (fInFormat == NULL)
196 return kERROR;
197
198 *fLog << (telrequired?inf:all) << "Open file: '" << name << "'" << endl;
199
200 if (telrequired)
201 {
202 // Check if thsi is a single telescope from a CORSIKA telescope array
203 // This is indicated by the -telXXX attached to the filename
204 TObjArray *res = TPRegexp("(.*/)?cer([0-9]+)-tel([0-9]+)").MatchS(name, "i");
205 if (res->GetLast()==3)
206 {
207 const TString match = res->At(0)->GetName();
208 const TString path = res->At(1)->GetName();
209 const UInt_t runid = atoi(res->At(2)->GetName());
210 const UInt_t telid = atoi(res->At(3)->GetName());
211
212 *fLog << inf2 << "Run " << runid << " detected to be for telescope #" << telid << endl;
213
214 if (telid>=fNumTelescopes)
215 {
216 *fLog << err << "Position of telescope " << telid << " not defined." << endl;
217 return kERROR;
218 }
219
220 *fLog << inf << "Telescope #" << telid << " [X/Y/Z (R)]: ";
221 *fLog << fTelescopeX[telid] << "/";
222 *fLog << fTelescopeY[telid] << "/";
223 *fLog << fTelescopeZ[telid] << " (R=" << fTelescopeR[telid] << ")" << endl;
224
225 fNumTelescope = telid+1;
226 }
227 delete res;
228 }
229
230 if (fDisplay)
231 {
232 // Show the number of the last event after
233 // which we now open a new file
234 TString txt = GetFileName();
235 txt += " @ ";
236 txt += GetNumExecutions()-1;
237 fDisplay->SetStatusLine2(txt);
238 }
239
240 fNumFile++;
241
242 if (!fParList)
243 return kTRUE;
244
245 //
246 // Search for MTaskList
247 //
248 MTask *tlist = (MTask*)fParList->FindObject("MTaskList");
249 if (!tlist)
250 {
251 *fLog << err << dbginf << "MTaskList not found... abort." << endl;
252 return kERROR;
253 }
254
255 //
256 // A new file has been opened and new headers have been read.
257 // --> ReInit tasklist
258 //
259 return tlist->ReInit(fParList) ? kTRUE : kERROR;
260}
261
262// --------------------------------------------------------------------------
263//
264// Return file name of current file.
265//
266TString MCorsikaRead::GetFullFileName() const
267{
268 const TObject *file = fFileNames->At(fNumFile-1);
269 return file ? file->GetName() : "";
270}
271
272// --------------------------------------------------------------------------
273//
274// Restart with the first file
275//
276Bool_t MCorsikaRead::Rewind()
277{
278 fNumFile=0;
279 fNumEvents=0;
280 return OpenNextFile()==kTRUE;
281}
282
283// --------------------------------------------------------------------------
284//
285// Calculates the total number of events from all input files and stores
286// the number in fNumTotalEvents.
287//
288Bool_t MCorsikaRead::CalcNumTotalEvents()
289{
290 fNumTotalEvents = 0;
291
292 bool allocated = false;
293 MCorsikaRunHeader *runheader = fRunHeader;
294 if (!runheader)
295 {
296 runheader = new MCorsikaRunHeader;
297 allocated = true;
298 }
299
300 Bool_t rc = kTRUE;
301 while (1)
302 {
303 switch (OpenNextFile(kFALSE, kFALSE))
304 {
305 case kFALSE:
306 break;
307 case kERROR:
308 rc = kFALSE;
309 break;
310 case kTRUE:
311
312 // read run header(1200), telescope position(1201) and
313 // first event header(1202)
314 Bool_t status = kTRUE;
315 Int_t blockType, blockVersion, blockIdentifier, blockLength;
316 while (status &&
317 fInFormat->NextBlock(fReadState/* == 4*/, blockType, blockVersion,
318 blockIdentifier, blockLength))
319 {
320 if (blockType == 1200)
321 status = runheader->ReadEvt(fInFormat, blockLength / sizeof(Float_t));
322
323 else if(blockType == 1201)
324 status = ReadTelescopePosition();
325
326 else if (blockType == 1202)
327 {
328 vector<Float_t> buffer(blockLength / sizeof(Float_t));
329 status = fInFormat->Read(buffer.data(), blockLength);
330 status = runheader->ReadEventHeader(buffer.data());
331 break;
332 }
333 else
334 fInFormat->Seek(blockLength);
335 }
336
337 if (status != kTRUE)
338 return status;
339
340 if (!fInFormat->SeekEvtEnd())
341 {
342 *fLog << (fForceMode?warn:err) << "ERROR - RUNE section not found in file." << endl;
343 if (allocated)
344 delete runheader;
345 if (!fForceMode)
346 return fForceMode ? kTRUE : kFALSE;
347 }
348
349 if (!runheader->ReadEvtEnd(fInFormat, kTRUE))
350 {
351 *fLog << (fForceMode?warn:err) << "ERROR - Reading RUNE section failed." << endl;
352 if (allocated)
353 delete runheader;
354 if (!fForceMode)
355 return kFALSE;
356 }
357
358 fNumTotalEvents += runheader->GetNumEvents()*runheader->GetNumReuse()*
359 (fTelescopeIdx<0 && fNumTelescope==0 ? fNumTelescopes : 1);
360 continue;
361 }
362 break;
363 }
364
365 if (allocated)
366 delete runheader;
367 return rc;
368}
369
370// --------------------------------------------------------------------------
371//
372// Reads the the position of all telescopes in one array
373//
374Int_t MCorsikaRead::ReadTelescopePosition()
375{
376 if (!fInFormat->Read(&fNumTelescopes, 4)) return kERROR;
377
378 if (fTelescopeIdx>=Int_t(fNumTelescopes))
379 {
380 *fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx <<
381 " exceeds number of telescopes " << fNumTelescopes << " in file." << endl;
382 return kERROR;
383 }
384
385 fTelescopeX.Set(fNumTelescopes);
386 fTelescopeY.Set(fNumTelescopes);
387 fTelescopeZ.Set(fNumTelescopes);
388 fTelescopeR.Set(fNumTelescopes);
389
390 if (!fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes)) return kERROR;
391 if (!fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes)) return kERROR;
392 if (!fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes)) return kERROR;
393 if (!fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes)) return kERROR;
394
395 *fLog << all;
396 *fLog << "Number of telescopes: " << fNumTelescopes << " (using ";
397 if (fTelescopeIdx>=0)
398 *fLog << "telescope " << fTelescopeIdx;
399 else
400 *fLog << "all telescopes";
401 *fLog << ")" << endl;
402
403 const Int_t lo = fTelescopeIdx<0 ? 0 : fTelescopeIdx;
404 const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1;
405
406 for (int i=lo; i<up; i++)
407 {
408 *fLog << all << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
409 *fLog << setw(7) << fTelescopeX[i] << "/";
410 *fLog << setw(7) << fTelescopeY[i] << "/";
411 *fLog << setw(7) << fTelescopeZ[i] << " (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
412 }
413
414 fNumTelescope = 0;
415
416 return kTRUE;
417}
418
419// --------------------------------------------------------------------------
420//
421// Reads the header of the next block. Depending on the current fReadState
422// and the block-type of the new header the methods decides wether
423// to continue the reading (kCONTINEU), to exit the Process() method (kTRUE)
424// or to try to read a new file (kFALSE).
425// Return codes:
426// - kFALSE : end of file was detected and no new header was read
427// - kTRU: A EventEnd was already found (fReadStatus == 3) and
428// the current header is not the RunEnd
429// - kCONTINUE: all other cases.
430Int_t MCorsikaRead::ReadNextBlockHeader()
431{
432 if (fInFormat->NextBlock(fReadState/* == 4*/, fBlockType, fBlockVersion,
433 fBlockIdentifier, fBlockLength) == kFALSE)
434 // end of file
435 return kFALSE;
436/*
437 gLog << "Next fBlock: type=" << fBlockType << " version=" << fBlockVersion;
438 gLog << " identifier=" << fBlockIdentifier << " length=" << fBlockLength;
439 gLog << " readState= " << fReadState << endl;
440*/
441 if (fReadState == 3 && fBlockType != 1210)
442 // fReadState == 3 means we have read the event end
443 // most of the time we can now save our data. BUT:
444 // bBlockType != 1210 means the next block is not RUNE,
445 // we want to read the RUNE block, before we
446 // save everything (last events and RUNE)
447 return kTRUE;
448
449 if (fBlockType == 1204 && fReadState != 2)
450 // next is a new set of telescope arrays, we store the previous ones
451 // but not if this is the first one (fReadState != 2)
452 return kTRUE;
453
454 return kCONTINUE;
455
456}
457
458// --------------------------------------------------------------------------
459//
460// The PreProcess of this task checks for the following containers in the
461// list:
462// MCorsikaRunHeader <output> if not found it is created
463// MCorsikaEvtHeader <output> if not found it is created
464// MCorsikaEvtData <output> if not found it is created
465// MCorsikaCrateArray <output> if not found it is created
466// MCorsikaEvtTime <output> if not found it is created (MTime)
467//
468// If all containers are found or created the run header is read from the
469// binary file and printed. If the Magic-Number (file identification)
470// doesn't match we stop the eventloop.
471//
472// Now the EvtHeader and EvtData containers are initialized.
473//
474Int_t MCorsikaRead::PreProcess(MParList *pList)
475{
476 //
477 // open the input stream
478 // first of all check if opening the file in the constructor was
479 // successfull
480 //
481 fParList = 0;
482
483 //
484 // check if all necessary containers exist in the Parameter list.
485 // if not create one and add them to the list
486 //
487 fRunHeader = (MCorsikaRunHeader*)pList->FindCreateObj("MCorsikaRunHeader");
488 if (!fRunHeader)
489 return kFALSE;
490
491 fEvtHeader = (MCorsikaEvtHeader*)pList->FindCreateObj("MCorsikaEvtHeader");
492 if (!fEvtHeader)
493 return kFALSE;
494
495 fEvent = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent");
496 if (!fEvent)
497 return kFALSE;
498
499 *fLog << inf << "Calculating number of total events..." << flush;
500 if (!CalcNumTotalEvents())
501 return kFALSE;
502 *fLog << inf << " " << fNumTotalEvents << " found." << endl;
503
504 fNumFile=0;
505 fNumEvents=0;
506
507 fParList = pList;
508
509 return kTRUE;
510}
511
512// --------------------------------------------------------------------------
513//
514// The Process reads one event from the binary file:
515// - The event header is read
516// - the run header is read (only once per file)
517// - the raw data information of one event is read
518//
519Int_t MCorsikaRead::Process()
520{
521 while (1) // loop for multiple input files
522 {
523 if (fInFormat)
524 {
525
526 while (1) // loop per input block
527 {
528
529 if (fReadState == 4)
530 {
531 fTopBlockLength -= fBlockLength + 12;
532 if (fTopBlockLength <= 0)
533 // all data of a top block are read, go back to normal state
534 fReadState = 15; // not used
535 }
536
537
538 Int_t status = kTRUE;
539 switch (fBlockType)
540 {
541 case 1200: // the run header
542 status = fRunHeader->ReadEvt(fInFormat, fBlockLength / sizeof(Float_t));
543 fRunHeader->SetTelescopeNumber(fNumTelescope);
544 fReadState = 1; // RUNH is read
545 break;
546
547 case 1201: // telescope positions
548 status = ReadTelescopePosition();
549 break;
550
551 case 1202: // the event header
552 {
553 vector<Float_t> buffer(fBlockLength / sizeof(Float_t));
554 if (!fInFormat->Read(buffer.data(), fBlockLength))
555 return kFALSE;
556
557 if (fReadState == 1) // first event after RUN header
558 {
559 fRunHeader->ReadEventHeader(buffer.data());
560 fRunHeader->Print();
561 }
562
563 status = fEvtHeader->ReadEvt(buffer.data());
564 if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
565 {
566 *fLog << err << "ERROR - Requested array index " << fArrayIdx <<
567 " exceeds number of arrays " << fEvtHeader->GetTotReuse() <<
568 " in file." << endl;
569 return kERROR;
570 }
571 }
572 fReadState = 2;
573 break;
574
575 case 1203: // 16 bytes
576 fInFormat->Seek(fBlockLength);
577 break;
578
579 case 1212:
580 {
581 char *buf = new char[fBlockLength];
582 fInFormat->Read(buf, fBlockLength);
583 status = kTRUE;
584
585 char *ptr = buf;
586
587 unsigned int n = reinterpret_cast<unsigned int*>(ptr)[0];
588 ptr += 4;
589
590 *fLog << all << endl;
591
592 for (unsigned int i=0; i<n && ptr<buf+fBlockLength; i++)
593 {
594 unsigned short s = reinterpret_cast<unsigned short*>(ptr)[0];
595 ptr += 2;
596
597 *fLog << string(ptr, ptr+s) << '\n';
598 ptr += s;
599 }
600 *fLog << '\n' << endl;
601
602 delete [] buf;
603 }
604 break;
605
606
607 case 1204: // top level block for one array (only for eventio data)
608 if (fArrayIdx < 0 || fArrayIdx == fBlockIdentifier)
609 {
610 fReadState = 4;
611 fTopBlockLength = fBlockLength;
612 }
613 else
614 // skip this array of telescopes
615 fInFormat->Seek(fBlockLength);
616
617 break;
618
619
620 case 1205: // eventio data
621 {
622 Int_t telIdx = fBlockIdentifier % 1000;
623 if ((fBlockVersion == 0 || fBlockVersion == 1000) &&
624 (fTelescopeIdx < 0 || fTelescopeIdx == telIdx) )
625 {
626 status = fBlockVersion==0 ? fEvent->ReadEventIoEvt(fInFormat) : fEvent->ReadEventIoEvtCompact(fInFormat);
627
628 Int_t arrayIdx = fBlockIdentifier / 1000;
629 Float_t xArrOff, yArrOff;
630 fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
631 fEvtHeader->SetTelescopeOffset(arrayIdx,
632 xArrOff + fTelescopeY[telIdx],
633 yArrOff - fTelescopeX[telIdx] );
634 fEvent->AddXY(xArrOff + fTelescopeY[telIdx],
635 yArrOff - fTelescopeX[telIdx]);
636 fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
637 fRunHeader->GetWavelengthMax());
638 }
639 else
640 // skip this telescope
641 fInFormat->Seek(fBlockLength);
642 }
643 break;
644
645 case 1209: // the event end
646 status = fEvtHeader->ReadEvtEnd(fInFormat, fBlockLength / sizeof(Float_t));
647
648 if (fReadState == 10 || fReadState == 2)
649 {
650 // corsika events
651 fEvtHeader->ResetNumReuse();
652 fEvtHeader->InitXY();
653 if (fNumTelescope>0)
654 {
655 // Here, the impact has opposite sign -- I don't understand why
656 fEvtHeader->SetXY(+fTelescopeY[fNumTelescope-1]-fEvtHeader->GetX(),
657 -fTelescopeX[fNumTelescope-1]-fEvtHeader->GetY());
658 }
659 fBlockType = 1109; // save corsika events
660 continue;
661 }
662
663 fReadState = 3; // event closed, run still open
664 break;
665
666
667 case 1210: // the run end
668 status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
669 fNumEvents += fRunHeader->GetNumEvents();
670 //fRunHeader->SetReadyToSave();
671 fReadState = 5; // back to starting point
672 // this may be the last block in the file.
673 // We exit to write the header into the file
674 fBlockLength = 0;
675 fBlockType = 0; // not available type
676 return kTRUE;
677 break;
678
679
680 case 1105: // event block of raw format
681 if (fReadState == 2 || fReadState == 10)
682 {
683 Int_t oldSize = fRawEvemtBuffer.size();
684 fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
685 status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength);
686 fReadState = 10;
687 }
688 else
689 fInFormat->Seek(fBlockLength);
690 break;
691
692
693 case 1109: // save corsika events
694 fEvtHeader->InitXY();
695 if (fNumTelescope>0)
696 {
697 // Here, the impact has opposite sign -- I don't understand why
698 fEvtHeader->SetXY(+fTelescopeY[fNumTelescope-1]-fEvtHeader->GetX(),
699 -fTelescopeX[fNumTelescope-1]-fEvtHeader->GetY());
700 }
701 status = fEvent->ReadCorsikaEvt(fRawEvemtBuffer,
702 fBlockLength == MCorsikaFormat::kBlockLengthRaw/21 - 4 ? 7 : 8,
703 fEvtHeader->GetNumReuse()+1);
704
705 // Simulate wavelength for all bunches with a wavelength == 0
706 fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
707 fRunHeader->GetWavelengthMax());
708
709 fEvtHeader->IncNumReuse();
710
711 if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse())
712 {
713 // this was the last reuse. Set fBlockType to EVTE to save
714 // it the next time.
715 fRawEvemtBuffer.resize(0);
716
717 fReadState = 3;
718 fBlockType = 1209;
719 }
720 else
721 // save this reuse
722 return status;
723
724 break;
725
726 default:
727 // unknown block, we skip it
728 fReadState = 15;
729 fInFormat->Seek(fBlockLength);
730
731 }
732
733 if (status != kTRUE)
734 // there was an error while data were read
735 return status;
736
737 Int_t headerStatus = ReadNextBlockHeader();
738 if (headerStatus == kFALSE)
739 // end of file
740 break;
741 if (headerStatus == kTRUE)
742 // we shall quit this loop
743 return kTRUE;
744
745 // else continue
746 }
747
748 }
749
750 //
751 // If an event could not be read from file try to open new file
752 //
753 const Int_t rc = OpenNextFile();
754 if (rc!=kTRUE)
755 return rc;
756
757 if (ReadNextBlockHeader() == kFALSE)
758 return kFALSE;
759 }
760 return kTRUE;
761}
762
763
764// --------------------------------------------------------------------------
765//
766// Close the file. Check whether the number of read events differs from
767// the number the file should containe (MCorsikaRunHeader). Prints a warning
768// if it doesn't match.
769//
770Int_t MCorsikaRead::PostProcess()
771{
772
773 const UInt_t n = fNumTotalEvents;//fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1);
774
775 //
776 // Sanity check for the number of events
777 //
778 if (n==GetNumExecutions()-1 || GetNumExecutions()==0)
779 return kTRUE;
780
781 *fLog << warn << dec;
782 *fLog << "Warning - Number of read events (" << GetNumExecutions()-1;
783 *fLog << ") doesn't match number in RUNE section(s) (" << n << ")";
784 //*fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
785
786 return kTRUE;
787}
788
789Int_t MCorsikaRead::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
790{
791 Bool_t rc = kFALSE;
792 if (IsEnvDefined(env, prefix, "TelescopeX", print) || IsEnvDefined(env, prefix, "TelescopeY", print))
793 {
794 rc = kTRUE;
795
796 fTelescopeX.Set(1);
797 fTelescopeY.Set(1);
798 fNumTelescope = 1;
799
800 if (IsEnvDefined(env, prefix, "TelescopeX", print))
801 fTelescopeX[0] = GetEnvValue(env, prefix, "TelescopeX", 0.);
802 if (IsEnvDefined(env, prefix, "TelescopeY", print))
803 fTelescopeY[0] = GetEnvValue(env, prefix, "TelescopeY", 0.);
804
805 *fLog << all << "Telescope #" << setw(4) << fNumTelescope << " [X/Y]: ";
806 *fLog << setw(7) << fTelescopeX[0] << "/";
807 *fLog << setw(7) << fTelescopeY[0] << endl;
808 }
809
810 // Read telescope positions from corsika input card
811 if (IsEnvDefined(env, prefix, "CorsikaInputCard", print))
812 {
813 rc = kTRUE;
814 TString fname = GetEnvValue(env, prefix, "CorsikaInputCard", "");
815
816 gSystem->ExpandPathName(fname);
817
818 // Is file accessible
819 if (gSystem->AccessPathName(fname, kFileExists))
820 return kERROR;
821
822 *fLog << inf2 << "Reading Telecope positions from " << fname << ":" << endl;
823
824 ifstream fin(fname);
825
826
827 TString buf;
828 while (1)
829 {
830 buf.ReadLine(fin);
831 if (!fin)
832 break;
833
834 buf = buf.Strip(TString::kBoth);
835 if (buf(0, 9)!="TELESCOPE")
836 continue;
837
838 float x, y, z, r;
839 unsigned int idx;
840 if (sscanf(buf.Data()+9, "%f %f %f %f %ud", &x, &y, &z, &r, &idx)!=5)
841 continue;
842
843 if (idx>=fNumTelescopes)
844 {
845 fNumTelescopes = idx+1;
846 fTelescopeX.Set(fNumTelescopes);
847 fTelescopeY.Set(fNumTelescopes);
848 fTelescopeZ.Set(fNumTelescopes);
849 fTelescopeR.Set(fNumTelescopes);
850 }
851
852 fTelescopeX[idx] = x;
853 fTelescopeY[idx] = y;
854 fTelescopeZ[idx] = z;
855 fTelescopeR[idx] = r;
856 }
857
858 for (UInt_t i=0; i<fNumTelescopes; i++)
859 {
860 *fLog << all << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
861 *fLog << setw(7) << fTelescopeX[i] << "/";
862 *fLog << setw(7) << fTelescopeY[i] << "/";
863 *fLog << setw(7) << fTelescopeZ[i] << " (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
864 }
865 }
866
867 return rc;
868}
869
Note: See TracBrowser for help on using the repository browser.