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