Context Navigation

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

Visit:

Last change on this file since 19797 was 19696, checked in by tbretz, 5 years ago
Fixed the caluclation of num events in case of single telescope files belonging to multi telescope setups.
File size: 26.7 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
60	ClassImp(MCorsikaRead);
61
62	using namespace std;
63
64
65	/* ----------- please don't delete and don't care about (Thomas) ------------
66	#define kBUFSZ 64 //1024102464
67	#include <iomanip.h>
68	class bifstream : public istream, public streambuf
69	{
70	private:
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	}
88	public:
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	//
101	MCorsikaRead::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	//
120	MCorsikaRead::~MCorsikaRead()
121	{
122	delete fFileNames;
123	if (fInFormat)
124	delete fInFormat;
125	}
126
127	/*
128	Byte_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	//
158	Int_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	//
169	Int_t MCorsikaRead::OpenNextFile(Bool_t print)
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 << inf << "Open file: '" << name << "'" << endl;
199
200	// Check if thsi is a single telescope from a CORSIKA telescope array
201	// This is indicated by the -telXXX attached to the filename
202	TObjArray res = TPRegexp("(./)?cer([0-9]+)-tel([0-9]+)").MatchS(name, "i");
203	if (res->GetLast()==3)
204	{
205	const TString match = res->At(0)->GetName();
206	const TString path = res->At(1)->GetName();
207	const UInt_t runid = atoi(res->At(2)->GetName());
208	const UInt_t telid = atoi(res->At(3)->GetName());
209
210	*fLog << inf2 << "Run " << runid << " detected to be for telescope #" << telid << endl;
211
212	if (telid>=fNumTelescopes)
213	{
214	*fLog << err << "Position of telescope " << telid << " not defined." << endl;
215	return kERROR;
216	}
217
218	*fLog << inf << "Telescope #" << telid << " [X/Y/Z (R)]: ";
219	*fLog << fTelescopeX[telid] << "/";
220	*fLog << fTelescopeY[telid] << "/";
221	*fLog << fTelescopeZ[telid] << " (R=" << fTelescopeR[telid] << ")" << endl;
222
223	fNumTelescope = telid+1;
224	}
225	delete res;
226
227
228	if (fDisplay)
229	{
230	// Show the number of the last event after
231	// which we now open a new file
232	TString txt = GetFileName();
233	txt += " @ ";
234	txt += GetNumExecutions()-1;
235	fDisplay->SetStatusLine2(txt);
236	}
237
238	fNumFile++;
239
240	if (!fParList)
241	return kTRUE;
242
243	//
244	// Search for MTaskList
245	//
246	MTask tlist = (MTask)fParList->FindObject("MTaskList");
247	if (!tlist)
248	{
249	*fLog << err << dbginf << "MTaskList not found... abort." << endl;
250	return kERROR;
251	}
252
253	//
254	// A new file has been opened and new headers have been read.
255	// --> ReInit tasklist
256	//
257	return tlist->ReInit(fParList) ? kTRUE : kERROR;
258	}
259
260	// --------------------------------------------------------------------------
261	//
262	// Return file name of current file.
263	//
264	TString MCorsikaRead::GetFullFileName() const
265	{
266	const TObject *file = fFileNames->At(fNumFile-1);
267	return file ? file->GetName() : "";
268	}
269
270	// --------------------------------------------------------------------------
271	//
272	// Restart with the first file
273	//
274	Bool_t MCorsikaRead::Rewind()
275	{
276	fNumFile=0;
277	fNumEvents=0;
278	return OpenNextFile()==kTRUE;
279	}
280
281	// --------------------------------------------------------------------------
282	//
283	// Calculates the total number of events from all input files and stores
284	// the number in fNumTotalEvents.
285	//
286	Bool_t MCorsikaRead::CalcNumTotalEvents()
287	{
288	fNumTotalEvents = 0;
289
290	Bool_t rc = kTRUE;
291
292	while (1)
293	{
294	switch (OpenNextFile(kFALSE))
295	{
296	case kFALSE:
297	break;
298	case kERROR:
299	rc = kFALSE;
300	break;
301	case kTRUE:
302
303	// read run header(1200), telescope position(1201) and
304	// first event header(1202)
305	Bool_t status = kTRUE;
306	Int_t blockType, blockVersion, blockIdentifier, blockLength;
307	while (status &&
308	fInFormat->NextBlock(fReadState/* == 4*/, blockType, blockVersion,
309	blockIdentifier, blockLength))
310	{
311	if (blockType == 1200)
312	status = fRunHeader->ReadEvt(fInFormat, blockLength / sizeof(Float_t));
313
314	else if(blockType == 1201)
315	status = ReadTelescopePosition();
316
317	else if (blockType == 1202)
318	{
319	vector<Float_t> buffer(blockLength / sizeof(Float_t));
320	status = fInFormat->Read(buffer.data(), blockLength);
321	status = fRunHeader->ReadEventHeader(buffer.data());
322	break;
323	}
324	else
325	fInFormat->Seek(blockLength);
326	}
327
328	if (status != kTRUE)
329	return status;
330
331	if (!fInFormat->SeekEvtEnd())
332	{
333	*fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl;
334	if (!fForceMode)
335	return fForceMode ? kTRUE : kFALSE;
336	}
337
338	if (!fRunHeader->ReadEvtEnd(fInFormat, kTRUE))
339	{
340	*fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl;
341	if (!fForceMode)
342	return kFALSE;
343	}
344
345	fNumTotalEvents += fRunHeader->GetNumEvents()fRunHeader->GetNumReuse()
346	(fTelescopeIdx<0 && fNumTelescope==0 ? fNumTelescopes : 1);
347	continue;
348	}
349	break;
350	}
351
352	return rc;
353	}
354
355	// --------------------------------------------------------------------------
356	//
357	// Reads the the position of all telescopes in one array
358	//
359	Int_t MCorsikaRead::ReadTelescopePosition()
360	{
361	if (!fInFormat->Read(&fNumTelescopes, 4)) return kERROR;
362
363	if (fTelescopeIdx>=Int_t(fNumTelescopes))
364	{
365	*fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx <<
366	" exceeds number of telescopes " << fNumTelescopes << " in file." << endl;
367	return kERROR;
368	}
369
370	fTelescopeX.Set(fNumTelescopes);
371	fTelescopeY.Set(fNumTelescopes);
372	fTelescopeZ.Set(fNumTelescopes);
373	fTelescopeR.Set(fNumTelescopes);
374
375	if (!fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes)) return kERROR;
376	if (!fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes)) return kERROR;
377	if (!fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes)) return kERROR;
378	if (!fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes)) return kERROR;
379
380	*fLog << all;
381	*fLog << "Number of telescopes: " << fNumTelescopes << " (using ";
382	if (fTelescopeIdx>=0)
383	*fLog << "telescope " << fTelescopeIdx;
384	else
385	*fLog << "all telescopes";
386	*fLog << ")" << endl;
387
388	const Int_t lo = fTelescopeIdx<0 ? 0 : fTelescopeIdx;
389	const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1;
390
391	for (int i=lo; i<up; i++)
392	{
393	*fLog << all << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
394	*fLog << setw(7) << fTelescopeX[i] << "/";
395	*fLog << setw(7) << fTelescopeY[i] << "/";
396	*fLog << setw(7) << fTelescopeZ[i] << " (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
397	}
398
399	fNumTelescope = 0;
400
401	return kTRUE;
402	}
403
404	// --------------------------------------------------------------------------
405	//
406	// Reads the header of the next block. Depending on the current fReadState
407	// and the block-type of the new header the methods decides wether
408	// to continue the reading (kCONTINEU), to exit the Process() method (kTRUE)
409	// or to try to read a new file (kFALSE).
410	// Return codes:
411	// - kFALSE : end of file was detected and no new header was read
412	// - kTRU: A EventEnd was already found (fReadStatus == 3) and
413	// the current header is not the RunEnd
414	// - kCONTINUE: all other cases.
415	Int_t MCorsikaRead::ReadNextBlockHeader()
416	{
417	if (fInFormat->NextBlock(fReadState/* == 4*/, fBlockType, fBlockVersion,
418	fBlockIdentifier, fBlockLength) == kFALSE)
419	// end of file
420	return kFALSE;
421	/*
422	gLog << "Next fBlock: type=" << fBlockType << " version=" << fBlockVersion;
423	gLog << " identifier=" << fBlockIdentifier << " length=" << fBlockLength;
424	gLog << " readState= " << fReadState << endl;
425	*/
426	if (fReadState == 3 && fBlockType != 1210)
427	// fReadState == 3 means we have read the event end
428	// most of the time we can now save our data. BUT:
429	// bBlockType != 1210 means the next block is not RUNE,
430	// we want to read the RUNE block, before we
431	// save everything (last events and RUNE)
432	return kTRUE;
433
434	if (fBlockType == 1204 && fReadState != 2)
435	// next is a new set of telescope arrays, we store the previous ones
436	// but not if this is the first one (fReadState != 2)
437	return kTRUE;
438
439	return kCONTINUE;
440
441	}
442
443	// --------------------------------------------------------------------------
444	//
445	// The PreProcess of this task checks for the following containers in the
446	// list:
447	// MCorsikaRunHeader <output> if not found it is created
448	// MCorsikaEvtHeader <output> if not found it is created
449	// MCorsikaEvtData <output> if not found it is created
450	// MCorsikaCrateArray <output> if not found it is created
451	// MCorsikaEvtTime <output> if not found it is created (MTime)
452	//
453	// If all containers are found or created the run header is read from the
454	// binary file and printed. If the Magic-Number (file identification)
455	// doesn't match we stop the eventloop.
456	//
457	// Now the EvtHeader and EvtData containers are initialized.
458	//
459	Int_t MCorsikaRead::PreProcess(MParList *pList)
460	{
461	//
462	// open the input stream
463	// first of all check if opening the file in the constructor was
464	// successfull
465	//
466	fParList = 0;
467
468	//
469	// check if all necessary containers exist in the Parameter list.
470	// if not create one and add them to the list
471	//
472	fRunHeader = (MCorsikaRunHeader*)pList->FindCreateObj("MCorsikaRunHeader");
473	if (!fRunHeader)
474	return kFALSE;
475
476	fEvtHeader = (MCorsikaEvtHeader*)pList->FindCreateObj("MCorsikaEvtHeader");
477	if (!fEvtHeader)
478	return kFALSE;
479
480	fEvent = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent");
481	if (!fEvent)
482	return kFALSE;
483
484	*fLog << inf << "Calculating number of total events..." << flush;
485	if (!CalcNumTotalEvents())
486	return kFALSE;
487	*fLog << inf << " " << fNumTotalEvents << " found." << endl;
488
489	fNumFile=0;
490	fNumEvents=0;
491
492	fParList = pList;
493
494	return kTRUE;
495	}
496
497	// --------------------------------------------------------------------------
498	//
499	// The Process reads one event from the binary file:
500	// - The event header is read
501	// - the run header is read (only once per file)
502	// - the raw data information of one event is read
503	//
504	Int_t MCorsikaRead::Process()
505	{
506	while (1) // loop for multiple input files
507	{
508	if (fInFormat)
509	{
510
511	while (1) // loop per input block
512	{
513
514	if (fReadState == 4)
515	{
516	fTopBlockLength -= fBlockLength + 12;
517	if (fTopBlockLength <= 0)
518	// all data of a top block are read, go back to normal state
519	fReadState = 15; // not used
520	}
521
522
523	Int_t status = kTRUE;
524	switch (fBlockType)
525	{
526	case 1200: // the run header
527	status = fRunHeader->ReadEvt(fInFormat, fBlockLength / sizeof(Float_t));
528	fReadState = 1; // RUNH is read
529	break;
530
531	case 1201: // telescope positions
532	status = ReadTelescopePosition();
533	break;
534
535	case 1202: // the event header
536	{
537	vector<Float_t> buffer(fBlockLength / sizeof(Float_t));
538	if (!fInFormat->Read(buffer.data(), fBlockLength))
539	return kFALSE;
540
541	if (fReadState == 1) // first event after RUN header
542	{
543	fRunHeader->ReadEventHeader(buffer.data());
544	fRunHeader->Print();
545	}
546
547	status = fEvtHeader->ReadEvt(buffer.data());
548	if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
549	{
550	*fLog << err << "ERROR - Requested array index " << fArrayIdx <<
551	" exceeds number of arrays " << fEvtHeader->GetTotReuse() <<
552	" in file." << endl;
553	return kERROR;
554	}
555	}
556	fReadState = 2;
557	break;
558
559	case 1203: // 16 bytes
560	fInFormat->Seek(fBlockLength);
561	break;
562
563	case 1212:
564	{
565	char *buf = new char[fBlockLength];
566	fInFormat->Read(buf, fBlockLength);
567	status = kTRUE;
568
569	char *ptr = buf;
570
571	unsigned int n = reinterpret_cast<unsigned int*>(ptr)[0];
572	ptr += 4;
573
574	*fLog << all << endl;
575
576	for (unsigned int i=0; i<n && ptr<buf+fBlockLength; i++)
577	{
578	unsigned short s = reinterpret_cast<unsigned short*>(ptr)[0];
579	ptr += 2;
580
581	*fLog << string(ptr, ptr+s) << '\n';
582	ptr += s;
583	}
584	*fLog << '\n' << endl;
585
586	delete [] buf;
587	}
588	break;
589
590
591	case 1204: // top level block for one array (only for eventio data)
592	if (fArrayIdx < 0 \|\| fArrayIdx == fBlockIdentifier)
593	{
594	fReadState = 4;
595	fTopBlockLength = fBlockLength;
596	}
597	else
598	// skip this array of telescopes
599	fInFormat->Seek(fBlockLength);
600
601	break;
602
603
604	case 1205: // eventio data
605	{
606	Int_t telIdx = fBlockIdentifier % 1000;
607	if ((fBlockVersion == 0 \|\| fBlockVersion == 1000) &&
608	(fTelescopeIdx < 0 \|\| fTelescopeIdx == telIdx) )
609	{
610	status = fBlockVersion==0 ? fEvent->ReadEventIoEvt(fInFormat) : fEvent->ReadEventIoEvtCompact(fInFormat);
611
612	Int_t arrayIdx = fBlockIdentifier / 1000;
613	Float_t xArrOff, yArrOff;
614	fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
615	fEvtHeader->SetTelescopeOffset(arrayIdx,
616	xArrOff + fTelescopeY[telIdx],
617	yArrOff - fTelescopeX[telIdx] );
618	fEvent->AddXY(xArrOff + fTelescopeY[telIdx],
619	yArrOff - fTelescopeX[telIdx]);
620	fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
621	fRunHeader->GetWavelengthMax());
622	}
623	else
624	// skip this telescope
625	fInFormat->Seek(fBlockLength);
626	}
627	break;
628
629	case 1209: // the event end
630	status = fEvtHeader->ReadEvtEnd(fInFormat, fBlockLength / sizeof(Float_t));
631
632	if (fReadState == 10 \|\| fReadState == 2)
633	{
634	// corsika events
635	fEvtHeader->ResetNumReuse();
636	fEvtHeader->InitXY();
637	if (fNumTelescope>0)
638	{
639	// Here, the impact has opposite sign -- I don't understand why
640	fEvtHeader->SetXY(+fTelescopeY[fNumTelescope-1]-fEvtHeader->GetX(),
641	-fTelescopeX[fNumTelescope-1]-fEvtHeader->GetY());
642	}
643	fBlockType = 1109; // save corsika events
644	continue;
645	}
646
647	fReadState = 3; // event closed, run still open
648	break;
649
650
651	case 1210: // the run end
652	status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
653	fNumEvents += fRunHeader->GetNumEvents();
654	//fRunHeader->SetReadyToSave();
655	fReadState = 5; // back to starting point
656	// this may be the last block in the file.
657	// We exit to write the header into the file
658	fBlockLength = 0;
659	fBlockType = 0; // not available type
660	return kTRUE;
661	break;
662
663
664	case 1105: // event block of raw format
665	if (fReadState == 2 \|\| fReadState == 10)
666	{
667	Int_t oldSize = fRawEvemtBuffer.size();
668	fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
669	status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength);
670	fReadState = 10;
671	}
672	else
673	fInFormat->Seek(fBlockLength);
674	break;
675
676
677	case 1109: // save corsika events
678	fEvtHeader->InitXY();
679	if (fNumTelescope>0)
680	{
681	// Here, the impact has opposite sign -- I don't understand why
682	fEvtHeader->SetXY(+fTelescopeY[fNumTelescope-1]-fEvtHeader->GetX(),
683	-fTelescopeX[fNumTelescope-1]-fEvtHeader->GetY());
684	}
685	status = fEvent->ReadCorsikaEvt(fRawEvemtBuffer,
686	fBlockLength == MCorsikaFormat::kBlockLengthRaw/21 - 4 ? 7 : 8,
687	fEvtHeader->GetNumReuse()+1);
688
689	// Simulate wavelength for all bunches with a wavelength == 0
690	fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
691	fRunHeader->GetWavelengthMax());
692
693	fEvtHeader->IncNumReuse();
694
695	if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse())
696	{
697	// this was the last reuse. Set fBlockType to EVTE to save
698	// it the next time.
699	fRawEvemtBuffer.resize(0);
700
701	fReadState = 3;
702	fBlockType = 1209;
703	}
704	else
705	// save this reuse
706	return status;
707
708	break;
709
710	default:
711	// unknown block, we skip it
712	fReadState = 15;
713	fInFormat->Seek(fBlockLength);
714
715	}
716
717	if (status != kTRUE)
718	// there was an error while data were read
719	return status;
720
721	Int_t headerStatus = ReadNextBlockHeader();
722	if (headerStatus == kFALSE)
723	// end of file
724	break;
725	if (headerStatus == kTRUE)
726	// we shall quit this loop
727	return kTRUE;
728
729	// else continue
730	}
731
732	}
733
734	//
735	// If an event could not be read from file try to open new file
736	//
737	const Int_t rc = OpenNextFile();
738	if (rc!=kTRUE)
739	return rc;
740
741	if (ReadNextBlockHeader() == kFALSE)
742	return kFALSE;
743	}
744	return kTRUE;
745	}
746
747
748	// --------------------------------------------------------------------------
749	//
750	// Close the file. Check whether the number of read events differs from
751	// the number the file should containe (MCorsikaRunHeader). Prints a warning
752	// if it doesn't match.
753	//
754	Int_t MCorsikaRead::PostProcess()
755	{
756
757	const UInt_t n = fNumTotalEvents;//fNumEventsfRunHeader->GetNumReuse()(fTelescopeIdx<0?fNumTelescopes:1);
758
759	//
760	// Sanity check for the number of events
761	//
762	if (n==GetNumExecutions()-1 \|\| GetNumExecutions()==0)
763	return kTRUE;
764
765	*fLog << warn << dec;
766	*fLog << "Warning - Number of read events (" << GetNumExecutions()-1;
767	*fLog << ") doesn't match number in RUNE section(s) (" << n << ")";
768	//fLog << fRunHeader->GetNumReuse() << "" << fNumEvents << "=" << n << ")." << endl;
769
770	return kTRUE;
771	}
772
773	Int_t MCorsikaRead::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
774	{
775	Bool_t rc = kFALSE;
776	if (IsEnvDefined(env, prefix, "TelescopeX", print) \|\| IsEnvDefined(env, prefix, "TelescopeY", print))
777	{
778	rc = kTRUE;
779
780	fTelescopeX.Set(1);
781	fTelescopeY.Set(1);
782	fNumTelescope = 1;
783
784	if (IsEnvDefined(env, prefix, "TelescopeX", print))
785	fTelescopeX[0] = GetEnvValue(env, prefix, "TelescopeX", 0.);
786	if (IsEnvDefined(env, prefix, "TelescopeY", print))
787	fTelescopeY[0] = GetEnvValue(env, prefix, "TelescopeY", 0.);
788
789	*fLog << all << "Telescope #" << setw(4) << fNumTelescope << " [X/Y]: ";
790	*fLog << setw(7) << fTelescopeX[0] << "/";
791	*fLog << setw(7) << fTelescopeY[0] << endl;
792	}
793
794	// Read telescope positions from corsika input card
795	if (IsEnvDefined(env, prefix, "CorsikaInputCard", print))
796	{
797	rc = kTRUE;
798	TString fname = GetEnvValue(env, prefix, "CorsikaInputCard", "");
799
800	gSystem->ExpandPathName(fname);
801
802	// Is file accessible
803	if (gSystem->AccessPathName(fname, kFileExists))
804	return kERROR;
805
806	*fLog << inf2 << "Reading Telecope positions from " << fname << ":" << endl;
807
808	TPRegexp regexp("\\s*TELESCOPE"
809	"\\s+([-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?)"
810	"\\s+([-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?)"
811	"\\s+([-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?)"
812	"\\s+([-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?)"
813	"\\s+([0-9]+)\\s*");
814
815	ifstream fin(fname);
816
817
818	TString buf;
819	while (1)
820	{
821	buf.ReadLine(fin);
822	if (!fin)
823	break;
824
825	const TObjArray *res = regexp.MatchS(buf);
826	if (res->GetLast()==9)
827	{
828	const UInt_t idx = atoi(res->At(9)->GetName());
829	if (idx>=fNumTelescopes)
830	{
831	fNumTelescopes = idx+1;
832	fTelescopeX.Set(fNumTelescopes);
833	fTelescopeY.Set(fNumTelescopes);
834	fTelescopeZ.Set(fNumTelescopes);
835	fTelescopeR.Set(fNumTelescopes);
836	}
837
838	fTelescopeX[idx] = atof(res->At(1)->GetName());
839	fTelescopeY[idx] = atof(res->At(3)->GetName());
840	fTelescopeZ[idx] = atof(res->At(5)->GetName());
841	fTelescopeR[idx] = atof(res->At(7)->GetName());
842	}
843	delete res;
844	}
845
846	for (UInt_t i=0; i<fNumTelescopes; i++)
847	{
848	*fLog << all << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
849	*fLog << setw(7) << fTelescopeX[i] << "/";
850	*fLog << setw(7) << fTelescopeY[i] << "/";
851	*fLog << setw(7) << fTelescopeZ[i] << " (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
852	}
853	}
854
855	return rc;
856	}
857

Note: See TracBrowser for help on using the repository browser.

Download in other formats: