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