1 | #include "MonteCarlo.h"
|
---|
2 |
|
---|
3 | // --------------------------------------------------------------------------
|
---|
4 | // Default constructor. Initiates Pointers and Variables to NULL or 0
|
---|
5 | //
|
---|
6 | MonteCarlo::MonteCarlo()
|
---|
7 | {
|
---|
8 | InitVariables();
|
---|
9 |
|
---|
10 | return;
|
---|
11 | }
|
---|
12 |
|
---|
13 | MonteCarlo::MonteCarlo(
|
---|
14 | TString filename
|
---|
15 | )
|
---|
16 | {
|
---|
17 | InitVariables();
|
---|
18 |
|
---|
19 | mFileName = filename;
|
---|
20 |
|
---|
21 | // open root file return if it cannot be opend
|
---|
22 | if ( !OpenRootFile() ) return;
|
---|
23 |
|
---|
24 | LoadEventTree( "Events");
|
---|
25 | LoadHeaderTree( "RunHeaders");
|
---|
26 |
|
---|
27 | ReadRunHeader();
|
---|
28 |
|
---|
29 | return;
|
---|
30 | }
|
---|
31 |
|
---|
32 | MonteCarlo::MonteCarlo(
|
---|
33 | TString filename,
|
---|
34 | TString evtsTreeName
|
---|
35 | )
|
---|
36 | {
|
---|
37 | InitVariables();
|
---|
38 |
|
---|
39 | mFileName = filename;
|
---|
40 |
|
---|
41 | // open root file return if it cannot be opend
|
---|
42 | if ( !OpenRootFile() ) return;
|
---|
43 |
|
---|
44 | LoadEventTree( evtsTreeName);
|
---|
45 | LoadHeaderTree( "RunHeaders");
|
---|
46 |
|
---|
47 | ReadRunHeader();
|
---|
48 |
|
---|
49 | return;
|
---|
50 | }
|
---|
51 |
|
---|
52 | MonteCarlo::MonteCarlo(
|
---|
53 | TString filename,
|
---|
54 | TString evtsTreeName,
|
---|
55 | TString headerTreeName
|
---|
56 | )
|
---|
57 | {
|
---|
58 | InitVariables();
|
---|
59 |
|
---|
60 | mFileName = filename;
|
---|
61 |
|
---|
62 | // open root file return if it cannot be opend
|
---|
63 | if ( !OpenRootFile() ) return;
|
---|
64 |
|
---|
65 | LoadEventTree( evtsTreeName);
|
---|
66 | LoadHeaderTree( headerTreeName);
|
---|
67 |
|
---|
68 | ReadRunHeader();
|
---|
69 |
|
---|
70 | return;
|
---|
71 | }
|
---|
72 |
|
---|
73 | // --------------------------------------------------------------------------
|
---|
74 | // Destructor
|
---|
75 | //
|
---|
76 | MonteCarlo::~MonteCarlo(void)
|
---|
77 | {
|
---|
78 | delete[] mpPixel;
|
---|
79 |
|
---|
80 | CloseRootFile();
|
---|
81 |
|
---|
82 | delete mpRootFile;
|
---|
83 |
|
---|
84 | return;
|
---|
85 | }
|
---|
86 |
|
---|
87 | // --------------------------------------------------------------------------
|
---|
88 | // Initiates Pointers and Variables to NULL or 0
|
---|
89 | //
|
---|
90 | void
|
---|
91 | MonteCarlo::InitVariables()
|
---|
92 | {
|
---|
93 | // set variables' default values
|
---|
94 | mEventNumber = 0;
|
---|
95 | mNumberOfEvents = 2;
|
---|
96 | mVerbosityLvl = 3;
|
---|
97 | mFileName = "";
|
---|
98 |
|
---|
99 | // set default seperator in CSV
|
---|
100 | mSeparator = " ";
|
---|
101 |
|
---|
102 | // source file
|
---|
103 | mpRootFile = NULL;
|
---|
104 | mRootFileOpend = false;
|
---|
105 |
|
---|
106 | // Trees
|
---|
107 | mpEventTree = NULL;
|
---|
108 | mpHeaderTree = NULL;
|
---|
109 |
|
---|
110 |
|
---|
111 | //header data types
|
---|
112 | mpIntendedPulsePos = NULL;
|
---|
113 | mpMcRunHeader = NULL;
|
---|
114 | mpGeomCam = NULL;
|
---|
115 | mpRawRunHeader = NULL;
|
---|
116 | mpCorsikaRunHeader = NULL;
|
---|
117 |
|
---|
118 | //Evt data types
|
---|
119 | mpElectronicNoise = NULL;
|
---|
120 | mpRawEventData = NULL;
|
---|
121 | mpIncidentAngle = NULL;
|
---|
122 | mpMcEventMetaData = NULL;
|
---|
123 | mpRawEventHeader = NULL;
|
---|
124 | mpCorsikaEvtHeader = NULL;
|
---|
125 | //
|
---|
126 |
|
---|
127 | // containers
|
---|
128 | mpPixel = NULL;
|
---|
129 | // mpSamples = NULL;
|
---|
130 |
|
---|
131 | return;
|
---|
132 | }
|
---|
133 |
|
---|
134 | // ==========================================================================
|
---|
135 | // Setters
|
---|
136 | //
|
---|
137 |
|
---|
138 | void
|
---|
139 | MonteCarlo::SetVerbosityLevel(int verbLvl)
|
---|
140 | {
|
---|
141 | mVerbosityLvl = verbLvl;
|
---|
142 | return;
|
---|
143 | }
|
---|
144 |
|
---|
145 | void
|
---|
146 | MonteCarlo::SetSeparator(char sep)
|
---|
147 | {
|
---|
148 | mSeparator = sep;
|
---|
149 | return;
|
---|
150 | }
|
---|
151 |
|
---|
152 | // ==========================================================================
|
---|
153 | // Getters
|
---|
154 | //
|
---|
155 |
|
---|
156 | int
|
---|
157 | MonteCarlo::GetVerbosityLevel()
|
---|
158 | {
|
---|
159 | return mVerbosityLvl;
|
---|
160 | }
|
---|
161 |
|
---|
162 | char
|
---|
163 | MonteCarlo::GetSeparator()
|
---|
164 | {
|
---|
165 | return mSeparator;
|
---|
166 | }
|
---|
167 |
|
---|
168 | // ==========================================================================
|
---|
169 | // Root file handling
|
---|
170 | //
|
---|
171 |
|
---|
172 | int
|
---|
173 | MonteCarlo::OpenRootFile()
|
---|
174 | {
|
---|
175 | if (mVerbosityLvl > 0) cout << "...opening root file: " << mFileName << endl;
|
---|
176 |
|
---|
177 | mpRootFile = new TFile(mFileName, "READ");
|
---|
178 |
|
---|
179 | //check if root file could be opened
|
---|
180 | if (!mpRootFile->IsOpen())
|
---|
181 | {
|
---|
182 | cout << "ERROR - Could not find file " << mFileName << endl;
|
---|
183 | mRootFileOpend =false;
|
---|
184 | return 0;
|
---|
185 | }
|
---|
186 | mRootFileOpend = true;
|
---|
187 | return 1;
|
---|
188 | }
|
---|
189 |
|
---|
190 | void
|
---|
191 | MonteCarlo::CloseRootFile()
|
---|
192 | {
|
---|
193 | if (mVerbosityLvl > 0) cout << "...closing root file: " << mFileName << endl;
|
---|
194 |
|
---|
195 | // close root file
|
---|
196 | // If option == "R", all TProcessIDs referenced by this file are deleted.
|
---|
197 | mpRootFile->Close("R");
|
---|
198 |
|
---|
199 | mpRootFile=NULL;
|
---|
200 |
|
---|
201 | return;
|
---|
202 | }
|
---|
203 |
|
---|
204 | void
|
---|
205 | MonteCarlo::LoadHeaderTree(TString treeName)
|
---|
206 | {
|
---|
207 | if (mVerbosityLvl > 0) cout << "...loading run header tree" << endl;
|
---|
208 |
|
---|
209 | mpHeaderTree = (TTree*)mpRootFile->Get(treeName);
|
---|
210 |
|
---|
211 | //check if mpHeaderTree exists
|
---|
212 | if (mpHeaderTree->IsZombie())
|
---|
213 | {
|
---|
214 | cout << "...could not load tree" << endl;
|
---|
215 | return;
|
---|
216 | }
|
---|
217 |
|
---|
218 | // =======================================================================
|
---|
219 | //Set Adresses to Branches in RunHeader-Tree
|
---|
220 |
|
---|
221 | // -----------------------------------------------------------------------
|
---|
222 |
|
---|
223 | // MGeomCam
|
---|
224 | if ( mpHeaderTree->GetBranchStatus("MGeomCam.") )
|
---|
225 | {
|
---|
226 | if (mVerbosityLvl > 1) cout << " ...MGeomCam" << endl;
|
---|
227 | mpHeaderTree->SetBranchAddress("MGeomCam.", &mpGeomCam);
|
---|
228 | }
|
---|
229 | else cout << "...could not load branch: MGeomCam" << endl;
|
---|
230 |
|
---|
231 | // -----------------------------------------------------------------------
|
---|
232 |
|
---|
233 | // IntendedPulsePos
|
---|
234 | if ( mpHeaderTree->GetBranchStatus("IntendedPulsePos.") )
|
---|
235 | {
|
---|
236 | if (mVerbosityLvl > 1) cout << " ...IntendedPulsePos" << endl;
|
---|
237 | mpHeaderTree->SetBranchAddress("IntendedPulsePos.", &mpIntendedPulsePos);
|
---|
238 | }
|
---|
239 | else cout << "...could not load branch: IntendedPulsePos" << endl;
|
---|
240 |
|
---|
241 | // -----------------------------------------------------------------------
|
---|
242 |
|
---|
243 | // MMcRunHeader
|
---|
244 | if ( mpHeaderTree->GetBranchStatus("MMcRunHeader.") )
|
---|
245 | {
|
---|
246 | if (mVerbosityLvl > 1) cout << " ...MMcRunHeader" << endl;
|
---|
247 | mpHeaderTree->SetBranchAddress("MMcRunHeader.", &mpMcRunHeader);
|
---|
248 | }
|
---|
249 | else cout << "...could not load branch: MMcRunHeader" << endl;
|
---|
250 |
|
---|
251 | // -----------------------------------------------------------------------
|
---|
252 |
|
---|
253 | // ElectronicNoise
|
---|
254 | if ( mpHeaderTree->GetBranchStatus("ElectronicNoise.") )
|
---|
255 | {
|
---|
256 | if (mVerbosityLvl > 1) cout << " ...ElectronicNoise" << endl;
|
---|
257 | mpHeaderTree->SetBranchAddress("ElectronicNoise.", &mpElectronicNoise);
|
---|
258 | }
|
---|
259 | else cout << "...could not load branch: ElectronicNoise" << endl;
|
---|
260 |
|
---|
261 | // -----------------------------------------------------------------------
|
---|
262 |
|
---|
263 | // MRawRunHeader
|
---|
264 | if ( mpHeaderTree->GetBranchStatus("MRawRunHeader.") )
|
---|
265 | {
|
---|
266 | if (mVerbosityLvl > 1) cout << " ...MRawRunHeader" << endl;
|
---|
267 | mpHeaderTree->SetBranchAddress("MRawRunHeader.", &mpRawRunHeader);
|
---|
268 | }
|
---|
269 | else cout << "...could not load branch: MRawRunHeader" << endl;
|
---|
270 |
|
---|
271 | // -----------------------------------------------------------------------
|
---|
272 |
|
---|
273 | // MCorsikaRunHeader
|
---|
274 | if ( mpHeaderTree->GetBranchStatus("MCorsikaRunHeader.") )
|
---|
275 | {
|
---|
276 | if (mVerbosityLvl > 1) cout << " ...MCorsikaRunHeader" << endl;
|
---|
277 | mpHeaderTree->SetBranchAddress("MCorsikaRunHeader.",&mpCorsikaRunHeader);
|
---|
278 | }
|
---|
279 | else cout << "...could not load branch: MCorsikaRunHeader" << endl;
|
---|
280 |
|
---|
281 | // =======================================================================
|
---|
282 |
|
---|
283 | return;
|
---|
284 | }
|
---|
285 |
|
---|
286 |
|
---|
287 | void
|
---|
288 | MonteCarlo::ReadRunHeader()
|
---|
289 | {
|
---|
290 | // Read Values from RunHeader-Tree
|
---|
291 |
|
---|
292 | if (mVerbosityLvl > 0)
|
---|
293 | cout << "...reading run header " << mpHeaderTree << endl;
|
---|
294 |
|
---|
295 | //Get first and only entry
|
---|
296 | mpHeaderTree->GetEntry();
|
---|
297 |
|
---|
298 | // -----------------------------------------------------------------------
|
---|
299 |
|
---|
300 | //Get values from "MGeomCam" Branch
|
---|
301 | if ( mpGeomCam != NULL)
|
---|
302 | {
|
---|
303 | //Getter functions from "MGeomCamFACT.h"
|
---|
304 | mCamDist = mpGeomCam->GetCameraDist();
|
---|
305 | mNumberOfPixels = mpGeomCam->GetNumPixels();
|
---|
306 | mNumberOfSectors = mpGeomCam->GetNumSectors();
|
---|
307 | mNumberOfAreas = mpGeomCam->GetNumAreas();
|
---|
308 | }
|
---|
309 | else if (mVerbosityLvl > 2)
|
---|
310 | cout << "...could not read data from: MGeomCam" << endl;
|
---|
311 |
|
---|
312 | // -----------------------------------------------------------------------
|
---|
313 |
|
---|
314 | //Get values from "IntendedPulsePos" Branch
|
---|
315 | if ( mpIntendedPulsePos != NULL)
|
---|
316 | {
|
---|
317 | mIntendedPulsePos = mpIntendedPulsePos->GetVal();
|
---|
318 | }
|
---|
319 | else if (mVerbosityLvl > 2)
|
---|
320 | cout << "...could not read data from: IntendedPulsePos" << endl;
|
---|
321 |
|
---|
322 | // -----------------------------------------------------------------------
|
---|
323 |
|
---|
324 | //Get values from "MMcRunHeader" Branch from event Tree
|
---|
325 | if ( mpMcRunHeader != NULL)
|
---|
326 | {
|
---|
327 | //Getter functions from "MMcRunHeader.hxx"
|
---|
328 | mNumSimulatedShowers = mpMcRunHeader->GetNumSimulatedShowers();
|
---|
329 | }
|
---|
330 | else if (mVerbosityLvl > 2)
|
---|
331 | cout << "...could not read data from: MMcRunHeader" << endl;
|
---|
332 |
|
---|
333 | // -----------------------------------------------------------------------
|
---|
334 |
|
---|
335 | //Get number of Events from event Tree
|
---|
336 | if ( mpEventTree != NULL)
|
---|
337 | {
|
---|
338 | mNumberOfEntries = mpEventTree->GetEntries();
|
---|
339 | }
|
---|
340 | else if (mVerbosityLvl > 2)
|
---|
341 | cout << "...could not read number of Events from event Tree" << endl;
|
---|
342 |
|
---|
343 | if (mVerbosityLvl > 0)
|
---|
344 | cout << "Event Tree has " << mNumberOfEntries << " entries" << endl;
|
---|
345 |
|
---|
346 | // -----------------------------------------------------------------------
|
---|
347 |
|
---|
348 | //Get values from "MRawRunHeader" Branch
|
---|
349 | if ( mpRawRunHeader != NULL)
|
---|
350 | {
|
---|
351 | //Getter functions from "MRawRunHeader.h"
|
---|
352 | mNumberOfEvents = mpRawRunHeader->GetNumEvents();
|
---|
353 | mNumberOfEventsRead = mpRawRunHeader->GetNumEventsRead();
|
---|
354 | mSamplingFrequency = mpRawRunHeader->GetFreqSampling();
|
---|
355 | mSourceName = mpRawRunHeader->GetSourceName();
|
---|
356 | mFileNumber = mpRawRunHeader->GetFileNumber();
|
---|
357 | mNumberOfSamples = mpRawRunHeader->GetNumSamplesHiGain();
|
---|
358 | mNumberOfBytes = mpRawRunHeader->GetNumBytesPerSample();
|
---|
359 | mRunNumber = mpRawRunHeader->GetRunNumber();
|
---|
360 | mRunType = mpRawRunHeader->GetRunType();
|
---|
361 | }
|
---|
362 | else if (mVerbosityLvl > 2)
|
---|
363 | cout << "...could not read data from: MRawRunHeader" << endl;
|
---|
364 |
|
---|
365 | // -----------------------------------------------------------------------
|
---|
366 |
|
---|
367 | //Get values from "MCorsikaRunHeader" Branch
|
---|
368 | if ( mpCorsikaRunHeader != NULL)
|
---|
369 | {
|
---|
370 | //Getter functions from "MCorsikaRunHeader.h"
|
---|
371 | mSlopeSpectrum = mpCorsikaRunHeader->GetSlopeSpectrum();
|
---|
372 | mEnergyMin = mpCorsikaRunHeader->GetEnergyMin();
|
---|
373 | mEnergyMax = mpCorsikaRunHeader->GetEnergyMax();
|
---|
374 | mZdMin = mpCorsikaRunHeader->GetZdMin();
|
---|
375 | mZdMax = mpCorsikaRunHeader->GetZdMax();
|
---|
376 | mAzMin = mpCorsikaRunHeader->GetAzMin();
|
---|
377 | mAzMax = mpCorsikaRunHeader->GetAzMax();
|
---|
378 | }
|
---|
379 | else if (mVerbosityLvl > 2)
|
---|
380 | cout << "...could not read data from: MCorsikaRunHeader" << endl;
|
---|
381 |
|
---|
382 | // -----------------------------------------------------------------------
|
---|
383 |
|
---|
384 | // delete Pixel Array before you set refferences for a new one
|
---|
385 | // in case it is existing
|
---|
386 | if (mpPixel != NULL)
|
---|
387 | {
|
---|
388 | delete[] mpPixel;
|
---|
389 | }
|
---|
390 | mpPixel = new pixel_t[mNumberOfPixels];
|
---|
391 |
|
---|
392 | // -----------------------------------------------------------------------
|
---|
393 |
|
---|
394 | //Get Pedestal from "ElectronicNoise" Branch
|
---|
395 | if ( mpElectronicNoise != NULL)
|
---|
396 | {
|
---|
397 | if (mVerbosityLvl > 1) cout << " ...reading pedestal offsets" << endl;
|
---|
398 |
|
---|
399 | // Loop over all pixel: Read Pedestal Offset
|
---|
400 | for ( int i = 0; i < mNumberOfPixels; i++ )
|
---|
401 | {
|
---|
402 | // check if array entry exists
|
---|
403 | if ( &(mpElectronicNoise[0][i]) != NULL)
|
---|
404 | {
|
---|
405 | //tricky stuff!!!
|
---|
406 | // mpElectronicNoise is a MPedestalCam Array
|
---|
407 | // individual pixel pedestals are stored in a MPedestalPix array
|
---|
408 | // the [] operator is overloaded in MPedestalCam
|
---|
409 | // and returning a MPedestalPix at position [i]
|
---|
410 | // MPedestalPix hast a Getter called GetPedestal()
|
---|
411 | mpPixel[i].pedestal = mpElectronicNoise[0][i].GetPedestal();
|
---|
412 | }
|
---|
413 | else if (mVerbosityLvl > 2)
|
---|
414 | {
|
---|
415 | cout << " ...cannot read pedestal offset" << endl;
|
---|
416 | }
|
---|
417 | }
|
---|
418 | }
|
---|
419 | else if (mVerbosityLvl > 2)
|
---|
420 | cout << "...could not read data from: ElectronicNoise" << endl;
|
---|
421 |
|
---|
422 | return;
|
---|
423 | }
|
---|
424 |
|
---|
425 | void
|
---|
426 | MonteCarlo::LoadEventTree(TString treeName)
|
---|
427 | {
|
---|
428 | if (mVerbosityLvl > 0) cout << "...loading event tree" << endl;
|
---|
429 |
|
---|
430 | mpEventTree = (TTree*)mpRootFile->Get(treeName);
|
---|
431 |
|
---|
432 | if (mpEventTree->IsZombie())
|
---|
433 | {
|
---|
434 | cout << "...could not load tree" << endl;
|
---|
435 | return;
|
---|
436 | }
|
---|
437 |
|
---|
438 | // =======================================================================
|
---|
439 | //Set Adresses to Branches in Events-Tree
|
---|
440 |
|
---|
441 | if (mVerbosityLvl > 1) cout << "...SetBranchAddresses:" << endl;
|
---|
442 |
|
---|
443 | // -----------------------------------------------------------------------
|
---|
444 |
|
---|
445 | // MRawEvtData
|
---|
446 | if ( mpEventTree->GetBranchStatus("MRawEvtData.") != -1 )
|
---|
447 | {
|
---|
448 | if (mVerbosityLvl > 1) cout << " ...MRawEvtData" << endl;
|
---|
449 | mpEventTree->SetBranchAddress("MRawEvtData.", &mpRawEventData);
|
---|
450 | }
|
---|
451 | else cout << "...could not load branch: MRawEvtData" << endl;
|
---|
452 |
|
---|
453 | // -----------------------------------------------------------------------
|
---|
454 |
|
---|
455 | // IncidentAngle
|
---|
456 | if ( mpEventTree->GetBranchStatus("IncidentAngle.") )
|
---|
457 | {
|
---|
458 | //FIX ME: THIS VALUE IS NOT EXISTANT IN EVERY MC FILE
|
---|
459 |
|
---|
460 | if (mVerbosityLvl > 1) cout << " ...IncidentAngle" << endl;
|
---|
461 | mpEventTree->SetBranchAddress("IncidentAngle.", &mpIncidentAngle);
|
---|
462 | }
|
---|
463 | else cout << "...could not load branch: IncidentAngle" << endl;
|
---|
464 |
|
---|
465 | // -----------------------------------------------------------------------
|
---|
466 |
|
---|
467 | // MMcEvt
|
---|
468 | if ( mpEventTree->GetBranchStatus("MMcEvt.") )
|
---|
469 | {
|
---|
470 | if (mVerbosityLvl > 1) cout << " ...McEvt" << endl;
|
---|
471 | mpEventTree->SetBranchAddress("MMcEvt.", &mpMcEventMetaData);
|
---|
472 | }
|
---|
473 | else cout << "...could not load branch: McEvt" << endl;
|
---|
474 |
|
---|
475 | // -----------------------------------------------------------------------
|
---|
476 |
|
---|
477 | // MRawEvtHeader
|
---|
478 | if ( mpEventTree->GetBranchStatus("MRawEvtHeader.") )
|
---|
479 | {
|
---|
480 | if (mVerbosityLvl > 1) cout << " ...MRawEventHeader" << endl;
|
---|
481 | mpEventTree->SetBranchAddress("MRawEvtHeader.", &mpRawEventHeader);
|
---|
482 | }
|
---|
483 | else cout << "...could not load branch: MRawEvtHeader" << endl;
|
---|
484 |
|
---|
485 | // -----------------------------------------------------------------------
|
---|
486 |
|
---|
487 | // MCorsikaEvtHeader
|
---|
488 | if ( mpEventTree->GetBranchStatus("MCorsikaEvtHeader.") )
|
---|
489 | {
|
---|
490 | if (mVerbosityLvl > 1) cout << " ...MCorsikaEvtHeader" << endl;
|
---|
491 | mpEventTree->SetBranchAddress("MCorsikaEvtHeader.", &mpCorsikaEvtHeader);
|
---|
492 | }
|
---|
493 | else cout << "...could not load branch: MCorsikaEvtHeader" << endl;
|
---|
494 |
|
---|
495 | // =======================================================================
|
---|
496 |
|
---|
497 | return;
|
---|
498 | }
|
---|
499 |
|
---|
500 | void
|
---|
501 | MonteCarlo::ReadEventMetaData()
|
---|
502 | {
|
---|
503 | if (mVerbosityLvl > 1)
|
---|
504 | cout << "...reading event header" << endl;
|
---|
505 |
|
---|
506 | //Get values from "MGeomCamFACT" Branch
|
---|
507 | if ( mpIncidentAngle != NULL)
|
---|
508 | {
|
---|
509 | //Getter functions from "MGeomCamFACT.h"
|
---|
510 | mIncidentAngle = mpIncidentAngle->GetVal();
|
---|
511 | }
|
---|
512 | else if (mVerbosityLvl > 2)
|
---|
513 | cout << "...could not read data from: MGeomCamFACT" << endl;
|
---|
514 |
|
---|
515 | // -----------------------------------------------------------------------
|
---|
516 |
|
---|
517 | //Get values from "MMcEvt" Branch
|
---|
518 | if ( mpMcEventMetaData != NULL)
|
---|
519 | {
|
---|
520 | //The following Getter Functions can be found in MMcEvt.h
|
---|
521 | mCorsikaEventNumber = mpMcEventMetaData->GetEvtNumber();
|
---|
522 | mPhotElFromShower = mpMcEventMetaData->GetPhotElfromShower();
|
---|
523 | mPhotElinCamera = mpMcEventMetaData->GetPhotElinCamera();
|
---|
524 | //The following Getter Functions can be found in MMcEvtBasic.h
|
---|
525 | mPartId = mpMcEventMetaData->GetPartId();
|
---|
526 | mPartName = mpMcEventMetaData->GetParticleName(mPartId);
|
---|
527 | mPartSymbol = mpMcEventMetaData->GetParticleSymbol(mPartId);
|
---|
528 | mEnergy = mpMcEventMetaData->GetEnergy();
|
---|
529 | mImpact = mpMcEventMetaData->GetImpact();
|
---|
530 | mTelescopePhi = mpMcEventMetaData->GetTelescopePhi();
|
---|
531 | mTelescopeTheta = mpMcEventMetaData->GetTelescopeTheta();
|
---|
532 | mPhi = mpMcEventMetaData->GetParticlePhi();
|
---|
533 | mTheta = mpMcEventMetaData->GetParticleTheta();
|
---|
534 | }
|
---|
535 | else if (mVerbosityLvl > 2)
|
---|
536 | cout << "...could not read data from: MMcEvt" << endl;
|
---|
537 |
|
---|
538 | // -----------------------------------------------------------------------
|
---|
539 |
|
---|
540 | //Get values from "MRawEventHeader" Branch
|
---|
541 | if ( mpRawEventHeader != NULL)
|
---|
542 | {
|
---|
543 | //Getter functions from "MRawEventHeader.h"
|
---|
544 | mEventNumber = mpRawEventHeader->GetDAQEvtNumber();
|
---|
545 | mNumTriggerLvl1 = mpRawEventHeader->GetNumTrigLvl1();
|
---|
546 | mNumTriggerLvl2 = mpRawEventHeader->GetNumTrigLvl2();
|
---|
547 | }
|
---|
548 | else if (mVerbosityLvl > 2)
|
---|
549 | cout << "...could not read data from: MRawEventHeader" << endl;
|
---|
550 |
|
---|
551 | // -----------------------------------------------------------------------
|
---|
552 |
|
---|
553 | //Get values from "MCorsikaEvtHeader" Branch
|
---|
554 | if ( mpCorsikaEvtHeader != NULL)
|
---|
555 | {
|
---|
556 | //Getter functions from "MCorsikaEvtHeader.h"
|
---|
557 | mFirstInteractionHeight = mpCorsikaEvtHeader->GetFirstInteractionHeight();
|
---|
558 | mEvtReuse = mpCorsikaEvtHeader->GetNumReuse();
|
---|
559 | mMomentumX = mpCorsikaEvtHeader->GetMomentum().X();
|
---|
560 | mMomentumY = mpCorsikaEvtHeader->GetMomentum().Y();
|
---|
561 | mMomentumZ = mpCorsikaEvtHeader->GetMomentum().Z();
|
---|
562 | mZd = mpCorsikaEvtHeader->GetZd();
|
---|
563 | mAz = mpCorsikaEvtHeader->GetAz();
|
---|
564 | mX = mpCorsikaEvtHeader->GetX();
|
---|
565 | mY = mpCorsikaEvtHeader->GetY();
|
---|
566 | }
|
---|
567 | else if (mVerbosityLvl > 2)
|
---|
568 | cout << "...could not read data from: MCorsikaEvtHeader" << endl;
|
---|
569 |
|
---|
570 | // -----------------------------------------------------------------------
|
---|
571 |
|
---|
572 | // mWeightedNumPhotons = mpCorsikaEvtHeader->Get;
|
---|
573 | // no getter function, no idea how to extract information
|
---|
574 |
|
---|
575 | return;
|
---|
576 | }
|
---|
577 |
|
---|
578 | void
|
---|
579 | MonteCarlo::ReadEventRawData()
|
---|
580 | {
|
---|
581 | if (mVerbosityLvl > 1) cout << "...reading event raw data" << endl;
|
---|
582 |
|
---|
583 | // -----------------------------------------------------------------------
|
---|
584 |
|
---|
585 | // delete Pixel Array before you set refferences for a new one
|
---|
586 | // in case it is existing
|
---|
587 | if (mpPixel != NULL)
|
---|
588 | {
|
---|
589 | delete[] mpPixel;
|
---|
590 | }
|
---|
591 | mpPixel = new pixel_t[mNumberOfPixels];
|
---|
592 |
|
---|
593 | // -----------------------------------------------------------------------
|
---|
594 |
|
---|
595 | if ( mpRawEventData == NULL)
|
---|
596 | {
|
---|
597 | cout << "ERROR: cannot read event data" << endl;
|
---|
598 | return;
|
---|
599 | }
|
---|
600 |
|
---|
601 | // you have to set this before you can read information
|
---|
602 | // from a magic binary file
|
---|
603 | mpRawEventData->InitRead(mpRawRunHeader);
|
---|
604 |
|
---|
605 | int pix_first_sample;
|
---|
606 |
|
---|
607 | // -----------------------------------------------------------------------
|
---|
608 |
|
---|
609 | //array to contain alle samples of all pixel of a event
|
---|
610 | unsigned short* all_raw_data = NULL;
|
---|
611 |
|
---|
612 | if ( mpRawEventData != NULL)
|
---|
613 | {
|
---|
614 | // point raw data array to RawEventData Array in Event Tree
|
---|
615 | all_raw_data = (unsigned short*) mpRawEventData->GetSamples();
|
---|
616 | /*
|
---|
617 | // FADC samples (hi gain) of all pixels
|
---|
618 | // This is an array of Byte_t variables. The value of a FADC sample has a
|
---|
619 | // size of n=fNumBytesPerSample bytes. Therefore, a FADC sample value will
|
---|
620 | // occupy n consecutive elements of this array (little endian ordering, i.e,
|
---|
621 | // less significant bits (and bytes) go first.
|
---|
622 | // If m = GetNumHiGainSamples(), the n bytes corresponding to the value of the
|
---|
623 | // i-th FADC sample of the j-th pixel are stored in the n consecutive
|
---|
624 | // positions of this array, starting from fHiGainFadcSamples[j*n*m+i*n]
|
---|
625 | */
|
---|
626 | }
|
---|
627 | else cout << "...cannot read event raw data" << endl;
|
---|
628 |
|
---|
629 | // -----------------------------------------------------------------------
|
---|
630 |
|
---|
631 | if (mVerbosityLvl > 1)
|
---|
632 | cout << "...pixel progress: ";
|
---|
633 |
|
---|
634 | //Loop over all camera pixel and read eEvent raw data, pixel ID and pedestal offset
|
---|
635 | for ( int i = 0; i < mNumberOfPixels; i++ )
|
---|
636 | {
|
---|
637 | if (mVerbosityLvl > 1){
|
---|
638 | if ( !(i%(mNumberOfPixels/10) ))
|
---|
639 | cout << i/(mNumberOfPixels*1.0)*100 << "%.." ;
|
---|
640 | if ( i == mNumberOfPixels-1)
|
---|
641 | cout << "100% " ;
|
---|
642 | }
|
---|
643 |
|
---|
644 | // Read Pedestal Offset and store it in classes pixel_t array
|
---|
645 | if ( mpElectronicNoise != NULL)
|
---|
646 | {
|
---|
647 | mpPixel[i].pedestal = mpElectronicNoise[0][i].GetPedestal();
|
---|
648 | }
|
---|
649 | else cout << "...cannot read pedestal offset" << endl;
|
---|
650 |
|
---|
651 | // Read Pixel SoftId and store it in classes pixel_t array
|
---|
652 | if ( mpRawEventData != NULL)
|
---|
653 | {
|
---|
654 | mpPixel[i].SoftId = mpRawEventData->GetPixelIds()[i];
|
---|
655 | }
|
---|
656 | else cout << "...cannot read pixel id" << endl;
|
---|
657 |
|
---|
658 | pix_first_sample = i*mNumberOfSamples;
|
---|
659 |
|
---|
660 | // point beginning of class' pixel_t array to the address
|
---|
661 | // of pixel's first sample's adress in TTree
|
---|
662 | mpPixel[i].rawData = &(all_raw_data[pix_first_sample]);
|
---|
663 |
|
---|
664 | }
|
---|
665 |
|
---|
666 | if (mVerbosityLvl > 1)
|
---|
667 | cout << endl;
|
---|
668 |
|
---|
669 | return;
|
---|
670 | }
|
---|
671 |
|
---|
672 | void
|
---|
673 | MonteCarlo::ReadEvent(int Event)
|
---|
674 | {
|
---|
675 | if (mVerbosityLvl > 0) cout << endl
|
---|
676 | << "====================" << endl
|
---|
677 | << "...reading Event: " << Event << endl
|
---|
678 | << "====================" << endl;
|
---|
679 |
|
---|
680 |
|
---|
681 | //load certain event from tree
|
---|
682 | mpEventTree->GetEntry(Event);
|
---|
683 |
|
---|
684 | //Get Event header data from TTree
|
---|
685 | ReadEventMetaData();
|
---|
686 |
|
---|
687 | //Get Event raw data from TTree
|
---|
688 | ReadEventRawData();
|
---|
689 |
|
---|
690 | return;
|
---|
691 | }
|
---|
692 |
|
---|
693 | // ==========================================================================
|
---|
694 | // csv file handling
|
---|
695 | //
|
---|
696 |
|
---|
697 | void
|
---|
698 | MonteCarlo::OpenCsvFile(TString fileName)
|
---|
699 | {
|
---|
700 | mCsvFileName = fileName;
|
---|
701 |
|
---|
702 | if (mVerbosityLvl > 0) cout << "...opening csv file: " << mCsvFileName << endl;
|
---|
703 |
|
---|
704 | mCsv.open( mCsvFileName );
|
---|
705 |
|
---|
706 | return;
|
---|
707 | }
|
---|
708 |
|
---|
709 | void
|
---|
710 | MonteCarlo::CloseCsvFile()
|
---|
711 | {
|
---|
712 | if (mVerbosityLvl > 0) cout << "...closing csv file: " << mCsvFileName << endl;
|
---|
713 |
|
---|
714 | mCsv.close();
|
---|
715 |
|
---|
716 | return;
|
---|
717 | }
|
---|
718 |
|
---|
719 | void
|
---|
720 | MonteCarlo::WriteMc2Csv(TString filename)
|
---|
721 | {
|
---|
722 | if ( !mRootFileOpend ){
|
---|
723 | cout << "ERROR - no rootfile loaded, no data loaded, cannot write csv"
|
---|
724 | << endl;
|
---|
725 | return;
|
---|
726 | }
|
---|
727 |
|
---|
728 | cout << "...writing mc to csv: " << filename << endl;
|
---|
729 | cout << "...processing " << mNumberOfEntries << "Events" << endl;
|
---|
730 |
|
---|
731 | mCsvFileName = filename;
|
---|
732 | OpenCsvFile(mCsvFileName);
|
---|
733 |
|
---|
734 | WriteFileInfo2Csv();
|
---|
735 | WriteRunHeaderInfo2Csv();
|
---|
736 | WriteRunHeader2Csv();
|
---|
737 |
|
---|
738 | WritePedestalInfo2Csv();
|
---|
739 | WritePedestal2Csv();
|
---|
740 |
|
---|
741 | WriteEventHeaderInfo2Csv();
|
---|
742 | WriteEventDataInfo2Csv();
|
---|
743 |
|
---|
744 | // loop over all events from tree and write content to csv
|
---|
745 | // for (int evt = 0; evt < mNumberOfEvents; evt++)
|
---|
746 | for (int evt = 0; evt < mNumberOfEntries; evt++)
|
---|
747 | {
|
---|
748 | ReadEvent(evt);
|
---|
749 | WriteEventHeader2Csv();
|
---|
750 | WriteEventData2Csv();
|
---|
751 | }
|
---|
752 |
|
---|
753 | cout << endl << "...conversion done " << endl;
|
---|
754 |
|
---|
755 | CloseCsvFile();
|
---|
756 | return;
|
---|
757 | }
|
---|
758 |
|
---|
759 | void
|
---|
760 | MonteCarlo::WriteFileInfo2Csv()
|
---|
761 | {
|
---|
762 | if (mVerbosityLvl > 0) cout << "...writing file header to csv" << endl;
|
---|
763 |
|
---|
764 | mCsv << "### ==========================================================="
|
---|
765 | << endl;
|
---|
766 | mCsv << "### = FACT Monte Carlo ="
|
---|
767 | << endl;
|
---|
768 | mCsv << "### ==========================================================="
|
---|
769 | << endl;
|
---|
770 | mCsv << "### = FileInfo: "
|
---|
771 | << endl;
|
---|
772 | mCsv << "### = "
|
---|
773 | << endl;
|
---|
774 | mCsv << "### = Source Name: " << mSourceName << endl;
|
---|
775 | mCsv << "### = Number of Entries: " << mNumberOfEntries << endl;
|
---|
776 | mCsv << "### = Run Number: " << mRunNumber << endl;
|
---|
777 | mCsv << "### = Run Type : " << mRunType << endl;
|
---|
778 | mCsv << "### = File Number: " << mFileNumber << endl;
|
---|
779 | mCsv << "### ==========================================================="
|
---|
780 | << endl ;
|
---|
781 | mCsv << "###"
|
---|
782 | << endl ;
|
---|
783 |
|
---|
784 | return;
|
---|
785 | }
|
---|
786 |
|
---|
787 | void
|
---|
788 | MonteCarlo::WriteRunHeaderInfo2Csv()
|
---|
789 | {
|
---|
790 | if (mVerbosityLvl > 0) cout << "...writing run header names to csv" << endl;
|
---|
791 |
|
---|
792 | mCsv << "### [RunHeader]" << endl;
|
---|
793 |
|
---|
794 | mCsv << "# mNumberOfEntries" << mSeparator;
|
---|
795 | mCsv << "mIntendedPulsePos" << mSeparator;
|
---|
796 | // Csv << "mPedestalOffset" << mSeparator;
|
---|
797 | mCsv << "mNumSimulatedShowers" << mSeparator;
|
---|
798 | mCsv << "mNumberOfPixels" << mSeparator;
|
---|
799 | mCsv << "mNumberOfSamples" << mSeparator;
|
---|
800 | // mCsv << "mSampleSize" << mSeparator;
|
---|
801 | mCsv << "mCamDist" << mSeparator;
|
---|
802 | mCsv << "mSourceName" << mSeparator;
|
---|
803 | mCsv << "mSlopeSpectrum" << mSeparator;
|
---|
804 | mCsv << "mEnergyMin" << mSeparator;
|
---|
805 | mCsv << "mEnergyMax" << mSeparator;
|
---|
806 | mCsv << "mZdMin" << mSeparator;
|
---|
807 | mCsv << "mZdMax" << mSeparator;
|
---|
808 | mCsv << "mAzMin" << mSeparator;
|
---|
809 | mCsv << "mAzMax" << mSeparator;
|
---|
810 | mCsv << "mFileNumber" << mSeparator;
|
---|
811 | mCsv << "mRunnumber" << mSeparator;
|
---|
812 | mCsv << "mRunType" << endl;
|
---|
813 |
|
---|
814 | return;
|
---|
815 | }
|
---|
816 |
|
---|
817 | void
|
---|
818 | MonteCarlo::WriteRunHeader2Csv()
|
---|
819 | {
|
---|
820 | if (mVerbosityLvl > 0) cout << "...writing run header to csv" << endl;
|
---|
821 |
|
---|
822 | mCsv << mNumberOfEntries << mSeparator;
|
---|
823 | mCsv << mIntendedPulsePos << mSeparator;
|
---|
824 | // mCsv << mPedestalOffset << mSeparator;
|
---|
825 | mCsv << mNumSimulatedShowers << mSeparator;
|
---|
826 | mCsv << mNumberOfPixels << mSeparator;
|
---|
827 | mCsv << mNumberOfSamples << mSeparator;
|
---|
828 | // mCsv << mSampleSize << mSeparator;
|
---|
829 | mCsv << mCamDist << mSeparator;
|
---|
830 | mCsv << mSourceName << mSeparator;
|
---|
831 | mCsv << mSlopeSpectrum << mSeparator;
|
---|
832 | mCsv << mEnergyMin << mSeparator;
|
---|
833 | mCsv << mEnergyMax << mSeparator;
|
---|
834 | mCsv << mZdMin << mSeparator;
|
---|
835 | mCsv << mZdMax << mSeparator;
|
---|
836 | mCsv << mAzMin << mSeparator;
|
---|
837 | mCsv << mAzMax << mSeparator;
|
---|
838 | mCsv << mFileNumber << mSeparator;
|
---|
839 | mCsv << mRunNumber << mSeparator;
|
---|
840 | mCsv << mRunType << endl;
|
---|
841 |
|
---|
842 | return;
|
---|
843 | }
|
---|
844 |
|
---|
845 | void
|
---|
846 | MonteCarlo::WriteEventHeaderInfo2Csv()
|
---|
847 | {
|
---|
848 | if (mVerbosityLvl > 0) cout << "...writing event header names to csv" << endl;
|
---|
849 |
|
---|
850 | mCsv << "### [EventHeader]" << endl;
|
---|
851 |
|
---|
852 | mCsv << "# mEventNumber" << mSeparator;
|
---|
853 | mCsv << "mNumberOfBytes" << mSeparator;
|
---|
854 | mCsv << "mIncidentAngle" << mSeparator;
|
---|
855 | mCsv << "mPartId" << mSeparator;
|
---|
856 | mCsv << "mEnergy" << mSeparator;
|
---|
857 | mCsv << "mImpact" << mSeparator;
|
---|
858 | mCsv << "mTelescopePhi" << mSeparator;
|
---|
859 | mCsv << "mTelescopeTheta" << mSeparator;
|
---|
860 | mCsv << "mPhi" << mSeparator;
|
---|
861 | mCsv << "mTheta" << mSeparator;
|
---|
862 | mCsv << "mCorsikaEventNumber" << mSeparator;
|
---|
863 | mCsv << "mPhotElFromShower" << mSeparator;
|
---|
864 | mCsv << "mEvtReuse" << mSeparator;
|
---|
865 | mCsv << "mNumTriggerLvl1" << mSeparator;
|
---|
866 | mCsv << "mFirstInteractionHeight" << mSeparator;
|
---|
867 | mCsv << "mMomentumX" << mSeparator;
|
---|
868 | mCsv << "mMomentumY" << mSeparator;
|
---|
869 | mCsv << "mMomentumZ" << mSeparator;
|
---|
870 | mCsv << "mZd" << mSeparator;
|
---|
871 | mCsv << "mAz" << mSeparator;
|
---|
872 | mCsv << "mX" << mSeparator;
|
---|
873 | mCsv << "mY" << mSeparator;
|
---|
874 | // mCsv << "mWeightedNumPhotons" ;
|
---|
875 | mCsv << endl;
|
---|
876 |
|
---|
877 | return;
|
---|
878 | }
|
---|
879 |
|
---|
880 | void
|
---|
881 | MonteCarlo::WriteEventHeader2Csv()
|
---|
882 | {
|
---|
883 | if (mVerbosityLvl > 0) cout << "...writing event header to csv" << endl;
|
---|
884 |
|
---|
885 | mCsv << mEventNumber << mSeparator;
|
---|
886 | mCsv << mNumberOfBytes << mSeparator;
|
---|
887 | mCsv << mIncidentAngle << mSeparator;
|
---|
888 | mCsv << mPartId << mSeparator;
|
---|
889 | mCsv << mEnergy << mSeparator;
|
---|
890 | mCsv << mImpact << mSeparator;
|
---|
891 | mCsv << mTelescopePhi << mSeparator;
|
---|
892 | mCsv << mTelescopeTheta << mSeparator;
|
---|
893 | mCsv << mPhi << mSeparator;
|
---|
894 | mCsv << mTheta << mSeparator;
|
---|
895 | mCsv << mCorsikaEventNumber << mSeparator;
|
---|
896 | mCsv << mPhotElFromShower << mSeparator;
|
---|
897 | mCsv << mEvtReuse << mSeparator;
|
---|
898 | mCsv << mNumTriggerLvl1 << mSeparator;
|
---|
899 | mCsv << mFirstInteractionHeight << mSeparator;
|
---|
900 | mCsv << mMomentumX << mSeparator;
|
---|
901 | mCsv << mMomentumY << mSeparator;
|
---|
902 | mCsv << mMomentumZ << mSeparator;
|
---|
903 | mCsv << mZd << mSeparator;
|
---|
904 | mCsv << mAz << mSeparator;
|
---|
905 | mCsv << mX << mSeparator;
|
---|
906 | mCsv << mY << mSeparator;
|
---|
907 | // mCsv << mWeightedNumPhotons ;
|
---|
908 | mCsv << endl;
|
---|
909 |
|
---|
910 | return;
|
---|
911 | }
|
---|
912 |
|
---|
913 | void
|
---|
914 | MonteCarlo::WriteEventDataInfo2Csv()
|
---|
915 | {
|
---|
916 | if (mVerbosityLvl > 0) cout << "...writing event categories to csv" << endl;
|
---|
917 |
|
---|
918 | mCsv << "### [RawData]" << endl;
|
---|
919 | mCsv << "# mEventNumber" << mSeparator;
|
---|
920 | mCsv << "pixelID" << mSeparator;
|
---|
921 | // mCsv << "pixelPedestal" << mSeparator;
|
---|
922 |
|
---|
923 | // Loop over number of all samples in pixel
|
---|
924 | for (int i = 0; i < mNumberOfSamples; i++)
|
---|
925 | {
|
---|
926 | mCsv << "Raw_" << i << mSeparator;
|
---|
927 | }
|
---|
928 | mCsv << endl;
|
---|
929 |
|
---|
930 | return;
|
---|
931 | }
|
---|
932 |
|
---|
933 | void
|
---|
934 | MonteCarlo::WriteEventData2Csv()
|
---|
935 | {
|
---|
936 | if (mVerbosityLvl > 0) cout << "...writing event data to csv" << endl;
|
---|
937 |
|
---|
938 | // Loop over all pixels and write pixeldata to csv
|
---|
939 | for (int i = 0; i < mNumberOfPixels; i++)
|
---|
940 | {
|
---|
941 | WritePixelData2Csv(i);
|
---|
942 | }
|
---|
943 |
|
---|
944 | return;
|
---|
945 | }
|
---|
946 |
|
---|
947 | void
|
---|
948 | MonteCarlo::WritePixelData2Csv(int pixelID)
|
---|
949 | {
|
---|
950 | if (mVerbosityLvl > 3) cout << "...writing pixel data to csv" << endl;
|
---|
951 | mCsv << mEventNumber << mSeparator;
|
---|
952 | mCsv << mpPixel[pixelID].SoftId << mSeparator;
|
---|
953 | // mCsv << mpPixel[pixelID].pedestal << mSeparator;
|
---|
954 |
|
---|
955 | // Loop over number of all samples in pixel and write samples content to csv
|
---|
956 | for (int i = 0; i < mNumberOfSamples; i++)
|
---|
957 | {
|
---|
958 | mCsv << mpPixel[pixelID].rawData[i] << mSeparator;
|
---|
959 | }
|
---|
960 | mCsv << endl;
|
---|
961 |
|
---|
962 | return;
|
---|
963 | }
|
---|
964 |
|
---|
965 | void
|
---|
966 | MonteCarlo::WritePedestalInfo2Csv()
|
---|
967 | {
|
---|
968 | mCsv << "### [Pedestal]" << endl;
|
---|
969 | // mCsv << "# SoftID" ;
|
---|
970 | // mCsv << mSeparator ;
|
---|
971 | // mCsv << "Pedestal" << endl;
|
---|
972 | }
|
---|
973 |
|
---|
974 | void
|
---|
975 | MonteCarlo::WritePedestal2Csv()
|
---|
976 | {
|
---|
977 | if (mVerbosityLvl > 3) cout << "...writing pedestal info to csv" << endl;
|
---|
978 |
|
---|
979 | // Loop over all pixels and write pixel pedestal to csv
|
---|
980 | for (int pixelID = 0; pixelID < mNumberOfPixels; pixelID++)
|
---|
981 | {
|
---|
982 | // mCsv << mpPixel[pixelID].SoftId;
|
---|
983 | mCsv << mpPixel[pixelID].pedestal;
|
---|
984 | if (pixelID < mNumberOfPixels -1)
|
---|
985 | mCsv << mSeparator;
|
---|
986 | }
|
---|
987 | mCsv << endl;
|
---|
988 |
|
---|
989 | return;
|
---|
990 | }
|
---|
991 |
|
---|
992 |
|
---|
993 |
|
---|
994 |
|
---|
995 |
|
---|
996 |
|
---|
997 |
|
---|
998 |
|
---|