Context Navigation

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

Visit:

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

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

Download in other formats: