Context Navigation

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

Visit:

Last change on this file since 20031 was 19859, checked in by tbretz, 5 years ago
Fixed a compiler warning.
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, 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	//
266	TString 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	//
276	Bool_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	//
288	Bool_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	//
374	Int_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.
430	Int_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	//
474	Int_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	//
519	Int_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	fReadState = 1; // RUNH is read
544	break;
545
546	case 1201: // telescope positions
547	status = ReadTelescopePosition();
548	break;
549
550	case 1202: // the event header
551	{
552	vector<Float_t> buffer(fBlockLength / sizeof(Float_t));
553	if (!fInFormat->Read(buffer.data(), fBlockLength))
554	return kFALSE;
555
556	if (fReadState == 1) // first event after RUN header
557	{
558	fRunHeader->ReadEventHeader(buffer.data());
559	fRunHeader->Print();
560	}
561
562	status = fEvtHeader->ReadEvt(buffer.data());
563	if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
564	{
565	*fLog << err << "ERROR - Requested array index " << fArrayIdx <<
566	" exceeds number of arrays " << fEvtHeader->GetTotReuse() <<
567	" in file." << endl;
568	return kERROR;
569	}
570	}
571	fReadState = 2;
572	break;
573
574	case 1203: // 16 bytes
575	fInFormat->Seek(fBlockLength);
576	break;
577
578	case 1212:
579	{
580	char *buf = new char[fBlockLength];
581	fInFormat->Read(buf, fBlockLength);
582	status = kTRUE;
583
584	char *ptr = buf;
585
586	unsigned int n = reinterpret_cast<unsigned int*>(ptr)[0];
587	ptr += 4;
588
589	*fLog << all << endl;
590
591	for (unsigned int i=0; i<n && ptr<buf+fBlockLength; i++)
592	{
593	unsigned short s = reinterpret_cast<unsigned short*>(ptr)[0];
594	ptr += 2;
595
596	*fLog << string(ptr, ptr+s) << '\n';
597	ptr += s;
598	}
599	*fLog << '\n' << endl;
600
601	delete [] buf;
602	}
603	break;
604
605
606	case 1204: // top level block for one array (only for eventio data)
607	if (fArrayIdx < 0 \|\| fArrayIdx == fBlockIdentifier)
608	{
609	fReadState = 4;
610	fTopBlockLength = fBlockLength;
611	}
612	else
613	// skip this array of telescopes
614	fInFormat->Seek(fBlockLength);
615
616	break;
617
618
619	case 1205: // eventio data
620	{
621	Int_t telIdx = fBlockIdentifier % 1000;
622	if ((fBlockVersion == 0 \|\| fBlockVersion == 1000) &&
623	(fTelescopeIdx < 0 \|\| fTelescopeIdx == telIdx) )
624	{
625	status = fBlockVersion==0 ? fEvent->ReadEventIoEvt(fInFormat) : fEvent->ReadEventIoEvtCompact(fInFormat);
626
627	Int_t arrayIdx = fBlockIdentifier / 1000;
628	Float_t xArrOff, yArrOff;
629	fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
630	fEvtHeader->SetTelescopeOffset(arrayIdx,
631	xArrOff + fTelescopeY[telIdx],
632	yArrOff - fTelescopeX[telIdx] );
633	fEvent->AddXY(xArrOff + fTelescopeY[telIdx],
634	yArrOff - fTelescopeX[telIdx]);
635	fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
636	fRunHeader->GetWavelengthMax());
637	}
638	else
639	// skip this telescope
640	fInFormat->Seek(fBlockLength);
641	}
642	break;
643
644	case 1209: // the event end
645	status = fEvtHeader->ReadEvtEnd(fInFormat, fBlockLength / sizeof(Float_t));
646
647	if (fReadState == 10 \|\| fReadState == 2)
648	{
649	// corsika events
650	fEvtHeader->ResetNumReuse();
651	fEvtHeader->InitXY();
652	if (fNumTelescope>0)
653	{
654	// Here, the impact has opposite sign -- I don't understand why
655	fEvtHeader->SetXY(+fTelescopeY[fNumTelescope-1]-fEvtHeader->GetX(),
656	-fTelescopeX[fNumTelescope-1]-fEvtHeader->GetY());
657	}
658	fBlockType = 1109; // save corsika events
659	continue;
660	}
661
662	fReadState = 3; // event closed, run still open
663	break;
664
665
666	case 1210: // the run end
667	status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
668	fNumEvents += fRunHeader->GetNumEvents();
669	//fRunHeader->SetReadyToSave();
670	fReadState = 5; // back to starting point
671	// this may be the last block in the file.
672	// We exit to write the header into the file
673	fBlockLength = 0;
674	fBlockType = 0; // not available type
675	return kTRUE;
676	break;
677
678
679	case 1105: // event block of raw format
680	if (fReadState == 2 \|\| fReadState == 10)
681	{
682	Int_t oldSize = fRawEvemtBuffer.size();
683	fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
684	status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength);
685	fReadState = 10;
686	}
687	else
688	fInFormat->Seek(fBlockLength);
689	break;
690
691
692	case 1109: // save corsika events
693	fEvtHeader->InitXY();
694	if (fNumTelescope>0)
695	{
696	// Here, the impact has opposite sign -- I don't understand why
697	fEvtHeader->SetXY(+fTelescopeY[fNumTelescope-1]-fEvtHeader->GetX(),
698	-fTelescopeX[fNumTelescope-1]-fEvtHeader->GetY());
699	}
700	status = fEvent->ReadCorsikaEvt(fRawEvemtBuffer,
701	fBlockLength == MCorsikaFormat::kBlockLengthRaw/21 - 4 ? 7 : 8,
702	fEvtHeader->GetNumReuse()+1);
703
704	// Simulate wavelength for all bunches with a wavelength == 0
705	fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
706	fRunHeader->GetWavelengthMax());
707
708	fEvtHeader->IncNumReuse();
709
710	if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse())
711	{
712	// this was the last reuse. Set fBlockType to EVTE to save
713	// it the next time.
714	fRawEvemtBuffer.resize(0);
715
716	fReadState = 3;
717	fBlockType = 1209;
718	}
719	else
720	// save this reuse
721	return status;
722
723	break;
724
725	default:
726	// unknown block, we skip it
727	fReadState = 15;
728	fInFormat->Seek(fBlockLength);
729
730	}
731
732	if (status != kTRUE)
733	// there was an error while data were read
734	return status;
735
736	Int_t headerStatus = ReadNextBlockHeader();
737	if (headerStatus == kFALSE)
738	// end of file
739	break;
740	if (headerStatus == kTRUE)
741	// we shall quit this loop
742	return kTRUE;
743
744	// else continue
745	}
746
747	}
748
749	//
750	// If an event could not be read from file try to open new file
751	//
752	const Int_t rc = OpenNextFile();
753	if (rc!=kTRUE)
754	return rc;
755
756	if (ReadNextBlockHeader() == kFALSE)
757	return kFALSE;
758	}
759	return kTRUE;
760	}
761
762
763	// --------------------------------------------------------------------------
764	//
765	// Close the file. Check whether the number of read events differs from
766	// the number the file should containe (MCorsikaRunHeader). Prints a warning
767	// if it doesn't match.
768	//
769	Int_t MCorsikaRead::PostProcess()
770	{
771
772	const UInt_t n = fNumTotalEvents;//fNumEventsfRunHeader->GetNumReuse()(fTelescopeIdx<0?fNumTelescopes:1);
773
774	//
775	// Sanity check for the number of events
776	//
777	if (n==GetNumExecutions()-1 \|\| GetNumExecutions()==0)
778	return kTRUE;
779
780	*fLog << warn << dec;
781	*fLog << "Warning - Number of read events (" << GetNumExecutions()-1;
782	*fLog << ") doesn't match number in RUNE section(s) (" << n << ")";
783	//fLog << fRunHeader->GetNumReuse() << "" << fNumEvents << "=" << n << ")." << endl;
784
785	return kTRUE;
786	}
787
788	Int_t MCorsikaRead::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
789	{
790	Bool_t rc = kFALSE;
791	if (IsEnvDefined(env, prefix, "TelescopeX", print) \|\| IsEnvDefined(env, prefix, "TelescopeY", print))
792	{
793	rc = kTRUE;
794
795	fTelescopeX.Set(1);
796	fTelescopeY.Set(1);
797	fNumTelescope = 1;
798
799	if (IsEnvDefined(env, prefix, "TelescopeX", print))
800	fTelescopeX[0] = GetEnvValue(env, prefix, "TelescopeX", 0.);
801	if (IsEnvDefined(env, prefix, "TelescopeY", print))
802	fTelescopeY[0] = GetEnvValue(env, prefix, "TelescopeY", 0.);
803
804	*fLog << all << "Telescope #" << setw(4) << fNumTelescope << " [X/Y]: ";
805	*fLog << setw(7) << fTelescopeX[0] << "/";
806	*fLog << setw(7) << fTelescopeY[0] << endl;
807	}
808
809	// Read telescope positions from corsika input card
810	if (IsEnvDefined(env, prefix, "CorsikaInputCard", print))
811	{
812	rc = kTRUE;
813	TString fname = GetEnvValue(env, prefix, "CorsikaInputCard", "");
814
815	gSystem->ExpandPathName(fname);
816
817	// Is file accessible
818	if (gSystem->AccessPathName(fname, kFileExists))
819	return kERROR;
820
821	*fLog << inf2 << "Reading Telecope positions from " << fname << ":" << endl;
822
823	ifstream fin(fname);
824
825
826	TString buf;
827	while (1)
828	{
829	buf.ReadLine(fin);
830	if (!fin)
831	break;
832
833	buf = buf.Strip(TString::kBoth);
834	if (buf(0, 9)!="TELESCOPE")
835	continue;
836
837	float x, y, z, r;
838	unsigned int idx;
839	if (sscanf(buf.Data()+9, "%f %f %f %f %ud", &x, &y, &z, &r, &idx)!=5)
840	continue;
841
842	if (idx>=fNumTelescopes)
843	{
844	fNumTelescopes = idx+1;
845	fTelescopeX.Set(fNumTelescopes);
846	fTelescopeY.Set(fNumTelescopes);
847	fTelescopeZ.Set(fNumTelescopes);
848	fTelescopeR.Set(fNumTelescopes);
849	}
850
851	fTelescopeX[idx] = x;
852	fTelescopeY[idx] = y;
853	fTelescopeZ[idx] = z;
854	fTelescopeR[idx] = r;
855	}
856
857	for (UInt_t i=0; i<fNumTelescopes; i++)
858	{
859	*fLog << all << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
860	*fLog << setw(7) << fTelescopeX[i] << "/";
861	*fLog << setw(7) << fTelescopeY[i] << "/";
862	*fLog << setw(7) << fTelescopeZ[i] << " (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
863	}
864	}
865
866	return rc;
867	}
868

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

Download in other formats: