| 1 | #include "MonteCarlo.h"
|
|---|
| 2 |
|
|---|
| 3 | MonteCarlo::MonteCarlo()
|
|---|
| 4 | {
|
|---|
| 5 | InitVariables();
|
|---|
| 6 |
|
|---|
| 7 | return;
|
|---|
| 8 | }
|
|---|
| 9 |
|
|---|
| 10 | MonteCarlo::MonteCarlo(
|
|---|
| 11 | TString filename
|
|---|
| 12 | )
|
|---|
| 13 | {
|
|---|
| 14 | InitVariables();
|
|---|
| 15 |
|
|---|
| 16 | mFileName = filename;
|
|---|
| 17 |
|
|---|
| 18 | OpenRootFile();
|
|---|
| 19 |
|
|---|
| 20 | LoadEventTree( "Events");
|
|---|
| 21 | LoadHeaderTree( "RunHeaders");
|
|---|
| 22 |
|
|---|
| 23 | ReadRunHeader();
|
|---|
| 24 |
|
|---|
| 25 | return;
|
|---|
| 26 | }
|
|---|
| 27 |
|
|---|
| 28 | MonteCarlo::MonteCarlo(
|
|---|
| 29 | TString filename,
|
|---|
| 30 | TString evtsTreeName
|
|---|
| 31 | )
|
|---|
| 32 | {
|
|---|
| 33 | InitVariables();
|
|---|
| 34 |
|
|---|
| 35 | mFileName = filename;
|
|---|
| 36 |
|
|---|
| 37 | OpenRootFile();
|
|---|
| 38 |
|
|---|
| 39 | LoadEventTree( evtsTreeName);
|
|---|
| 40 | LoadHeaderTree( "RunHeaders");
|
|---|
| 41 |
|
|---|
| 42 | ReadRunHeader();
|
|---|
| 43 |
|
|---|
| 44 | return;
|
|---|
| 45 | }
|
|---|
| 46 |
|
|---|
| 47 | MonteCarlo::MonteCarlo(
|
|---|
| 48 | TString filename,
|
|---|
| 49 | TString evtsTreeName,
|
|---|
| 50 | TString headerTreeName
|
|---|
| 51 | )
|
|---|
| 52 | {
|
|---|
| 53 | InitVariables();
|
|---|
| 54 |
|
|---|
| 55 | mFileName = filename;
|
|---|
| 56 |
|
|---|
| 57 | OpenRootFile();
|
|---|
| 58 |
|
|---|
| 59 | LoadEventTree( evtsTreeName);
|
|---|
| 60 | LoadHeaderTree( headerTreeName);
|
|---|
| 61 |
|
|---|
| 62 | ReadRunHeader();
|
|---|
| 63 |
|
|---|
| 64 | return;
|
|---|
| 65 | }
|
|---|
| 66 |
|
|---|
| 67 | MonteCarlo::~MonteCarlo(void)
|
|---|
| 68 | {
|
|---|
| 69 | delete[] mpPixel;
|
|---|
| 70 |
|
|---|
| 71 | CloseRootFile();
|
|---|
| 72 |
|
|---|
| 73 | delete mpRootFile;
|
|---|
| 74 |
|
|---|
| 75 | return;
|
|---|
| 76 | }
|
|---|
| 77 |
|
|---|
| 78 | void
|
|---|
| 79 | MonteCarlo::InitVariables()
|
|---|
| 80 | {
|
|---|
| 81 | mFileName = "";
|
|---|
| 82 | mSeparator = " ";
|
|---|
| 83 |
|
|---|
| 84 | // source file
|
|---|
| 85 | mpRootFile = NULL;
|
|---|
| 86 |
|
|---|
| 87 | // Trees
|
|---|
| 88 | mpEventTree = NULL;
|
|---|
| 89 | mpHeaderTree = NULL;
|
|---|
| 90 |
|
|---|
| 91 |
|
|---|
| 92 | //header data types
|
|---|
| 93 | mpIntendedPulsePos = NULL;
|
|---|
| 94 | mpMcRunHeader = NULL;
|
|---|
| 95 | mpGeomCam = NULL;
|
|---|
| 96 | mpRawRunHeader = NULL;
|
|---|
| 97 | mpCorsikaRunHeader = NULL;
|
|---|
| 98 |
|
|---|
| 99 | //Evt data types
|
|---|
| 100 | mpElectronicNoise = NULL;
|
|---|
| 101 | mpRawEventData = NULL;
|
|---|
| 102 | mpIncidentAngle = NULL;
|
|---|
| 103 | mpMcEventMetaData = NULL;
|
|---|
| 104 | mpRawEventHeader = NULL;
|
|---|
| 105 | mpCorsikaEvtHeader = NULL;
|
|---|
| 106 | //
|
|---|
| 107 |
|
|---|
| 108 | // containers
|
|---|
| 109 | mpPixel = NULL;
|
|---|
| 110 | mpSamples = NULL;
|
|---|
| 111 | mEventNumber = 0;
|
|---|
| 112 | mNumberOfEvents = 2;
|
|---|
| 113 |
|
|---|
| 114 | mVerbosityLvl = 2;
|
|---|
| 115 |
|
|---|
| 116 | return;
|
|---|
| 117 | }
|
|---|
| 118 |
|
|---|
| 119 | void
|
|---|
| 120 | MonteCarlo::SetVerbosityLevel(int verbLvl)
|
|---|
| 121 | {
|
|---|
| 122 | mVerbosityLvl = verbLvl;
|
|---|
| 123 |
|
|---|
| 124 | return;
|
|---|
| 125 | }
|
|---|
| 126 |
|
|---|
| 127 | int
|
|---|
| 128 | MonteCarlo::GetVerbosityLevel()
|
|---|
| 129 | {
|
|---|
| 130 | return mVerbosityLvl;
|
|---|
| 131 | }
|
|---|
| 132 |
|
|---|
| 133 | void
|
|---|
| 134 | MonteCarlo::OpenRootFile()
|
|---|
| 135 | {
|
|---|
| 136 | if (mVerbosityLvl > 0) cout << "...opening root file: " << mFileName << endl;
|
|---|
| 137 |
|
|---|
| 138 | mpRootFile = new TFile(mFileName, "READ");
|
|---|
| 139 |
|
|---|
| 140 | if (!mpRootFile->IsOpen())
|
|---|
| 141 | {
|
|---|
| 142 | cout << "ERROR - Could not find file " << mFileName << endl;
|
|---|
| 143 | return ;
|
|---|
| 144 | }
|
|---|
| 145 |
|
|---|
| 146 | return;
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 | void
|
|---|
| 150 | MonteCarlo::CloseRootFile()
|
|---|
| 151 | {
|
|---|
| 152 | if (mVerbosityLvl > 0) cout << "...closing root file: " << mFileName << endl;
|
|---|
| 153 | mpRootFile->Close("R");
|
|---|
| 154 | mpRootFile=NULL;
|
|---|
| 155 |
|
|---|
| 156 | return;
|
|---|
| 157 | }
|
|---|
| 158 |
|
|---|
| 159 | void
|
|---|
| 160 | MonteCarlo::OpenCsvFile(TString fileName)
|
|---|
| 161 | {
|
|---|
| 162 | mCsvFileName = fileName;
|
|---|
| 163 |
|
|---|
| 164 | if (mVerbosityLvl > 0) cout << "...opening csv file: " << mCsvFileName << endl;
|
|---|
| 165 |
|
|---|
| 166 | mCsv.open( mCsvFileName );
|
|---|
| 167 |
|
|---|
| 168 | return;
|
|---|
| 169 | }
|
|---|
| 170 |
|
|---|
| 171 | void
|
|---|
| 172 | MonteCarlo::CloseCsvFile()
|
|---|
| 173 | {
|
|---|
| 174 | if (mVerbosityLvl > 0) cout << "...closing csv file: " << mCsvFileName << endl;
|
|---|
| 175 |
|
|---|
| 176 | mCsv.close();
|
|---|
| 177 |
|
|---|
| 178 | return;
|
|---|
| 179 | }
|
|---|
| 180 |
|
|---|
| 181 | void
|
|---|
| 182 | MonteCarlo::LoadHeaderTree(TString treeName)
|
|---|
| 183 | {
|
|---|
| 184 | if (mVerbosityLvl > 0) cout << "...loading run header tree" << endl;
|
|---|
| 185 |
|
|---|
| 186 | mpHeaderTree = (TTree*)mpRootFile->Get(treeName);
|
|---|
| 187 |
|
|---|
| 188 | if (mpHeaderTree->IsZombie())
|
|---|
| 189 | {
|
|---|
| 190 | cout << "...could not load tree" << endl;
|
|---|
| 191 | return;
|
|---|
| 192 | }
|
|---|
| 193 |
|
|---|
| 194 | //Set Adresses to Branches in RunHeader-Tree
|
|---|
| 195 | if ( mpHeaderTree->GetBranchStatus("MGeomCam.") )
|
|---|
| 196 | {
|
|---|
| 197 | if (mVerbosityLvl > 1) cout << " ...MGeomCam" << endl;
|
|---|
| 198 | mpHeaderTree->SetBranchAddress("MGeomCam.", &mpGeomCam);
|
|---|
| 199 | }
|
|---|
| 200 | if ( mpHeaderTree->GetBranchStatus("IntendedPulsePos.") )
|
|---|
| 201 | {
|
|---|
| 202 | if (mVerbosityLvl > 1) cout << " ...IntendedPulsePos" << endl;
|
|---|
| 203 | mpHeaderTree->SetBranchAddress("IntendedPulsePos.", &mpIntendedPulsePos);
|
|---|
| 204 | }
|
|---|
| 205 | if ( mpHeaderTree->GetBranchStatus("MMcRunHeader.") )
|
|---|
| 206 | {
|
|---|
| 207 | if (mVerbosityLvl > 1) cout << " ...MMcRunHeader" << endl;
|
|---|
| 208 | mpHeaderTree->SetBranchAddress("MMcRunHeader.", &mpMcRunHeader);
|
|---|
| 209 | }
|
|---|
| 210 | if ( mpHeaderTree->GetBranchStatus("ElectronicNoise.") )
|
|---|
| 211 | {
|
|---|
| 212 | if (mVerbosityLvl > 1) cout << " ...ElectronicNoise" << endl;
|
|---|
| 213 | mpHeaderTree->SetBranchAddress("ElectronicNoise.", &mpElectronicNoise);
|
|---|
| 214 | }
|
|---|
| 215 | if ( mpHeaderTree->GetBranchStatus("MRawRunHeader.") )
|
|---|
| 216 | {
|
|---|
| 217 | if (mVerbosityLvl > 1) cout << " ...MRawRunHeader" << endl;
|
|---|
| 218 | mpHeaderTree->SetBranchAddress("MRawRunHeader.", &mpRawRunHeader);
|
|---|
| 219 | }
|
|---|
| 220 | if ( mpHeaderTree->GetBranchStatus("MCorsikaRunHeader.") )
|
|---|
| 221 | {
|
|---|
| 222 | if (mVerbosityLvl > 1) cout << " ...MCorsikaRunHeader" << endl;
|
|---|
| 223 | mpHeaderTree->SetBranchAddress("MCorsikaRunHeader.",&mpCorsikaRunHeader);
|
|---|
| 224 | }
|
|---|
| 225 |
|
|---|
| 226 | return;
|
|---|
| 227 | }
|
|---|
| 228 |
|
|---|
| 229 | void
|
|---|
| 230 | MonteCarlo::ReadRunHeader()
|
|---|
| 231 | {
|
|---|
| 232 | if (mVerbosityLvl > 0)
|
|---|
| 233 | cout << "...reading run header " << mpHeaderTree << endl;
|
|---|
| 234 |
|
|---|
| 235 | //Read Values from RunHeader-Tree
|
|---|
| 236 |
|
|---|
| 237 | mpHeaderTree->GetEntry(); //This causes problems
|
|---|
| 238 |
|
|---|
| 239 |
|
|---|
| 240 | //Get values from "MGeomCam" Branch
|
|---|
| 241 | //Getter functions from "MGeomCamFACT.h"
|
|---|
| 242 | mCamDist = mpGeomCam->GetCameraDist();
|
|---|
| 243 | mNumberOfPixels = mpGeomCam->GetNumPixels();
|
|---|
| 244 | mNumberOfSectors = mpGeomCam->GetNumSectors();
|
|---|
| 245 | mNumberOfAreas = mpGeomCam->GetNumAreas();
|
|---|
| 246 |
|
|---|
| 247 | //Get values from "IntendedPulsePos" Branch
|
|---|
| 248 | mIntendedPulsePos = mpIntendedPulsePos->GetVal();
|
|---|
| 249 |
|
|---|
| 250 | //Get values from "MMcRunHeader" Branch from event Tree
|
|---|
| 251 | mNumSimulatedShowers = mpMcRunHeader->GetNumSimulatedShowers();
|
|---|
| 252 |
|
|---|
| 253 | //Get number of Events from event Tree
|
|---|
| 254 | mNumberOfEntries = mpEventTree->GetEntries();
|
|---|
| 255 | if (mVerbosityLvl > 0)
|
|---|
| 256 | cout << "Event Tree has " << mNumberOfEntries << " entries" << endl;
|
|---|
| 257 |
|
|---|
| 258 | //Get values from "MRawRunHeader" Branch
|
|---|
| 259 | mNumberOfEvents = mpRawRunHeader->GetNumEvents(); // couses problems
|
|---|
| 260 | mNumberOfEventsRead = mpRawRunHeader->GetNumEventsRead();
|
|---|
| 261 | mSamplingFrequency = mpRawRunHeader->GetFreqSampling();
|
|---|
| 262 | mSourceName = mpRawRunHeader->GetSourceName();
|
|---|
| 263 | mFileNumber = mpRawRunHeader->GetFileNumber();
|
|---|
| 264 | mNumberOfSamples = mpRawRunHeader->GetNumSamplesHiGain();
|
|---|
| 265 | mNumberOfBytes = mpRawRunHeader->GetNumBytesPerSample();
|
|---|
| 266 | mRunNumber = mpRawRunHeader->GetRunNumber();
|
|---|
| 267 | mRunType = mpRawRunHeader->GetRunType();
|
|---|
| 268 |
|
|---|
| 269 |
|
|---|
| 270 | //Get values from "MCorsikaRunHeader" Branch
|
|---|
| 271 | mSlopeSpectrum = mpCorsikaRunHeader->GetSlopeSpectrum();
|
|---|
| 272 | mEnergyMin = mpCorsikaRunHeader->GetEnergyMin();
|
|---|
| 273 | mEnergyMax = mpCorsikaRunHeader->GetEnergyMax();
|
|---|
| 274 | mZdMin = mpCorsikaRunHeader->GetZdMin();
|
|---|
| 275 | mZdMax = mpCorsikaRunHeader->GetZdMax();
|
|---|
| 276 | mAzMin = mpCorsikaRunHeader->GetAzMin();
|
|---|
| 277 | mAzMax = mpCorsikaRunHeader->GetAzMax();
|
|---|
| 278 |
|
|---|
| 279 | return;
|
|---|
| 280 | }
|
|---|
| 281 |
|
|---|
| 282 | void
|
|---|
| 283 | MonteCarlo::LoadEventTree(TString treeName)
|
|---|
| 284 | {
|
|---|
| 285 | if (mVerbosityLvl > 0) cout << "...loading event tree" << endl;
|
|---|
| 286 |
|
|---|
| 287 | mpEventTree = (TTree*)mpRootFile->Get(treeName);
|
|---|
| 288 |
|
|---|
| 289 | if (mpEventTree->IsZombie())
|
|---|
| 290 | {
|
|---|
| 291 | cout << "...could not load tree" << endl;
|
|---|
| 292 | return;
|
|---|
| 293 | }
|
|---|
| 294 |
|
|---|
| 295 |
|
|---|
| 296 | //Set Adresses to Branches in Events-Tree
|
|---|
| 297 | if (mVerbosityLvl > 1) cout << "...SetBranchAddresses:" << endl;
|
|---|
| 298 |
|
|---|
| 299 | if ( mpEventTree->GetBranchStatus("MRawEvtData.") != -1 )
|
|---|
| 300 | {
|
|---|
| 301 | if (mVerbosityLvl > 1) cout << " ...MRawEvtData" << endl;
|
|---|
| 302 | mpEventTree->SetBranchAddress("MRawEvtData.", &mpRawEventData);
|
|---|
| 303 | }
|
|---|
| 304 | if ( mpEventTree->GetBranchStatus("IncidentAngle.") )
|
|---|
| 305 | {
|
|---|
| 306 | //FIX ME: THIS VALUE IS NOT EXISTANT IN EVERY MC FILE
|
|---|
| 307 |
|
|---|
| 308 | if (mVerbosityLvl > 1) cout << " ...IncidentAngle" << endl;
|
|---|
| 309 | mpEventTree->SetBranchAddress("IncidentAngle.", &mpIncidentAngle);
|
|---|
| 310 | }
|
|---|
| 311 | if ( mpEventTree->GetBranchStatus("MMcEvt.") )
|
|---|
| 312 | {
|
|---|
| 313 | if (mVerbosityLvl > 1) cout << " ...McEvt" << endl;
|
|---|
| 314 | mpEventTree->SetBranchAddress("MMcEvt.", &mpMcEventMetaData);
|
|---|
| 315 | }
|
|---|
| 316 | if ( mpEventTree->GetBranchStatus("MRawEvtHeader.") )
|
|---|
| 317 | {
|
|---|
| 318 | if (mVerbosityLvl > 1) cout << " ...MRawEventHeader" << endl;
|
|---|
| 319 | mpEventTree->SetBranchAddress("MRawEvtHeader.", &mpRawEventHeader);
|
|---|
| 320 | }
|
|---|
| 321 | if ( mpEventTree->GetBranchStatus("MCorsikaEvtHeader.") )
|
|---|
| 322 | {
|
|---|
| 323 | if (mVerbosityLvl > 1) cout << " ...MCorsikaEvtHeader" << endl;
|
|---|
| 324 | mpEventTree->SetBranchAddress("MCorsikaEvtHeader.", &mpCorsikaEvtHeader);
|
|---|
| 325 | }
|
|---|
| 326 |
|
|---|
| 327 | return;
|
|---|
| 328 | }
|
|---|
| 329 |
|
|---|
| 330 | void
|
|---|
| 331 | MonteCarlo::ReadEventMetaData()
|
|---|
| 332 | {
|
|---|
| 333 | if (mVerbosityLvl > 1) cout << "...reading event header" << endl;
|
|---|
| 334 |
|
|---|
| 335 | //Get values from "MGeomCamFACT" Branch
|
|---|
| 336 | mIncidentAngle = mpIncidentAngle->GetVal();
|
|---|
| 337 |
|
|---|
| 338 | //Get values from "MMcEvt" Branch
|
|---|
| 339 | //The following Getter Functions can be found in MMcEvt.h
|
|---|
| 340 | mCorsikaEventNumber = mpMcEventMetaData->GetEvtNumber();
|
|---|
| 341 | mPhotElFromShower = mpMcEventMetaData->GetPhotElfromShower();
|
|---|
| 342 | mPhotElinCamera = mpMcEventMetaData->GetPhotElinCamera();
|
|---|
| 343 | //The following Getter Functions can be found in MMcEvtBasic.h
|
|---|
| 344 | mPartId = mpMcEventMetaData->GetPartId();
|
|---|
| 345 | mPartName = mpMcEventMetaData->GetParticleName(mPartId);
|
|---|
| 346 | mPartSymbol = mpMcEventMetaData->GetParticleSymbol(mPartId);
|
|---|
| 347 | mEnergy = mpMcEventMetaData->GetEnergy();
|
|---|
| 348 | mImpact = mpMcEventMetaData->GetImpact();
|
|---|
| 349 | mTelescopePhi = mpMcEventMetaData->GetTelescopePhi();
|
|---|
| 350 | mTelescopeTheta = mpMcEventMetaData->GetTelescopeTheta();
|
|---|
| 351 | mPhi = mpMcEventMetaData->GetParticlePhi();
|
|---|
| 352 | mTheta = mpMcEventMetaData->GetParticleTheta();
|
|---|
| 353 |
|
|---|
| 354 | //Get values from "MRawEventHeader" Branch
|
|---|
| 355 | mEventNumber = mpRawEventHeader->GetDAQEvtNumber();
|
|---|
| 356 | mNumTriggerLvl1 = mpRawEventHeader->GetNumTrigLvl1();
|
|---|
| 357 | mNumTriggerLvl2 = mpRawEventHeader->GetNumTrigLvl2();
|
|---|
| 358 |
|
|---|
| 359 | //Get values from "MCorsikaEvtHeader" Branch
|
|---|
| 360 | mFirstInteractionHeight = mpCorsikaEvtHeader->GetFirstInteractionHeight();
|
|---|
| 361 | mEvtReuse = mpCorsikaEvtHeader->GetNumReuse();
|
|---|
| 362 | mMomentumX = mpCorsikaEvtHeader->GetMomentum().X();
|
|---|
| 363 | mMomentumY = mpCorsikaEvtHeader->GetMomentum().Y();
|
|---|
| 364 | mMomentumZ = mpCorsikaEvtHeader->GetMomentum().Z();
|
|---|
| 365 | mZd = mpCorsikaEvtHeader->GetZd();
|
|---|
| 366 | mAz = mpCorsikaEvtHeader->GetAz();
|
|---|
| 367 | mX = mpCorsikaEvtHeader->GetX();
|
|---|
| 368 | mY = mpCorsikaEvtHeader->GetY();
|
|---|
| 369 |
|
|---|
| 370 | // mWeightedNumPhotons = mpCorsikaEvtHeader->Get;
|
|---|
| 371 | // no getter function, no idea how to extract information
|
|---|
| 372 |
|
|---|
| 373 | return;
|
|---|
| 374 | }
|
|---|
| 375 |
|
|---|
| 376 | void
|
|---|
| 377 | MonteCarlo::ReadEventRawData()
|
|---|
| 378 | {
|
|---|
| 379 | if (mVerbosityLvl > 1) cout << "...reading event raw data" << endl;
|
|---|
| 380 |
|
|---|
| 381 | // delete Pixel Array before you set refferences for a new one
|
|---|
| 382 | if (mpPixel != NULL)
|
|---|
| 383 | {
|
|---|
| 384 | delete[] mpPixel;
|
|---|
| 385 | }
|
|---|
| 386 | mpPixel = new pixel_t[mNumberOfPixels];
|
|---|
| 387 |
|
|---|
| 388 | mpRawEventData->InitRead(mpRawRunHeader);
|
|---|
| 389 |
|
|---|
| 390 | int pix_first_sample;
|
|---|
| 391 | // int pix_last_sample;
|
|---|
| 392 |
|
|---|
| 393 | unsigned short* all_raw_data = NULL;
|
|---|
| 394 |
|
|---|
| 395 | all_raw_data = (unsigned short*) mpRawEventData->GetSamples();
|
|---|
| 396 |
|
|---|
| 397 | if (mVerbosityLvl > 1)
|
|---|
| 398 | cout << "...pixel progress: ";
|
|---|
| 399 |
|
|---|
| 400 | for ( int i = 0; i < mNumberOfPixels; i++ )
|
|---|
| 401 | {
|
|---|
| 402 | if (mVerbosityLvl > 1){
|
|---|
| 403 | if ( !(i%(mNumberOfPixels/10) ))
|
|---|
| 404 | cout << i/(mNumberOfPixels*1.0)*100 << "%.." ;
|
|---|
| 405 | if ( i == mNumberOfPixels-1)
|
|---|
| 406 | cout << "100% " ;
|
|---|
| 407 | }
|
|---|
| 408 |
|
|---|
| 409 |
|
|---|
| 410 | mpPixel[i].pedestal = mpElectronicNoise[0][i].GetPedestal();
|
|---|
| 411 | mpPixel[i].SoftId = mpRawEventData->GetPixelIds()[i];
|
|---|
| 412 |
|
|---|
| 413 | pix_first_sample = i*mNumberOfSamples;
|
|---|
| 414 |
|
|---|
| 415 | // point beginning of pixel's rawdata array to address
|
|---|
| 416 | // of pixel's first sample's adress
|
|---|
| 417 |
|
|---|
| 418 | mpPixel[i].rawData = &(all_raw_data[pix_first_sample]);
|
|---|
| 419 |
|
|---|
| 420 | // for (int j = pix_first_sample; j < pix_last_sample; j++)
|
|---|
| 421 | // {
|
|---|
| 422 | // cout << mpPixel[i].rawData[j] << " ";
|
|---|
| 423 | // }
|
|---|
| 424 | // cout << endl;
|
|---|
| 425 |
|
|---|
| 426 | // // this could be improved by
|
|---|
| 427 | // DONE BUT NOT TESTED
|
|---|
| 428 | // for (int j = i*mSampleSize; j < (i+1)*mSampleSize; j++)
|
|---|
| 429 | // {
|
|---|
| 430 | // }
|
|---|
| 431 | }
|
|---|
| 432 |
|
|---|
| 433 | if (mVerbosityLvl > 1)
|
|---|
| 434 | cout << endl;
|
|---|
| 435 |
|
|---|
| 436 | return;
|
|---|
| 437 | }
|
|---|
| 438 |
|
|---|
| 439 | void
|
|---|
| 440 | MonteCarlo::ReadEvent(int Event)
|
|---|
| 441 | {
|
|---|
| 442 | if (mVerbosityLvl > 0) cout << endl
|
|---|
| 443 | << "====================" << endl
|
|---|
| 444 | << "...reading Event: " << Event << endl
|
|---|
| 445 | << "====================" << endl;
|
|---|
| 446 |
|
|---|
| 447 |
|
|---|
| 448 | //load certain event from tree
|
|---|
| 449 | mpEventTree->GetEntry(Event);
|
|---|
| 450 |
|
|---|
| 451 | ReadEventMetaData();
|
|---|
| 452 | ReadEventRawData();
|
|---|
| 453 |
|
|---|
| 454 | return;
|
|---|
| 455 | }
|
|---|
| 456 |
|
|---|
| 457 |
|
|---|
| 458 | void
|
|---|
| 459 | MonteCarlo::WriteMc2Csv(TString filename)
|
|---|
| 460 | {
|
|---|
| 461 | cout << "...writing mc to csv: " << filename << endl;
|
|---|
| 462 | cout << "...processing " << mNumberOfEntries << "Events" << endl;
|
|---|
| 463 |
|
|---|
| 464 | mCsvFileName = filename;
|
|---|
| 465 | OpenCsvFile(mCsvFileName);
|
|---|
| 466 |
|
|---|
| 467 | WriteFileHeader2Csv();
|
|---|
| 468 | WriteRunHeaderNames2Csv();
|
|---|
| 469 | WriteRunHeader2Csv();
|
|---|
| 470 | WriteEventHeaderNames2Csv();
|
|---|
| 471 | WriteEventDataNames2Csv();
|
|---|
| 472 |
|
|---|
| 473 | // for (int evt = 0; evt < mNumberOfEvents; evt++)
|
|---|
| 474 | for (int evt = 0; evt < mNumberOfEntries; evt++)
|
|---|
| 475 | {
|
|---|
| 476 | ReadEvent(evt);
|
|---|
| 477 | WriteEventHeader2Csv();
|
|---|
| 478 | WriteEventData2Csv();
|
|---|
| 479 | }
|
|---|
| 480 |
|
|---|
| 481 | cout << endl << "...conversion done " << endl;
|
|---|
| 482 |
|
|---|
| 483 | CloseCsvFile();
|
|---|
| 484 | return;
|
|---|
| 485 | }
|
|---|
| 486 |
|
|---|
| 487 |
|
|---|
| 488 | void
|
|---|
| 489 | MonteCarlo::WritePixelData2Csv(int pixelID)
|
|---|
| 490 | {
|
|---|
| 491 | if (mVerbosityLvl > 3) cout << "...writing pixel data to csv" << endl;
|
|---|
| 492 | mCsv << mEventNumber << mSeparator;
|
|---|
| 493 | mCsv << mpPixel[pixelID].SoftId << mSeparator;
|
|---|
| 494 | mCsv << mpPixel[pixelID].pedestal << mSeparator;
|
|---|
| 495 |
|
|---|
| 496 | for (int i = 0; i < mNumberOfSamples; i++)
|
|---|
| 497 | {
|
|---|
| 498 | mCsv << mpPixel[pixelID].rawData[i] << mSeparator;
|
|---|
| 499 | }
|
|---|
| 500 | mCsv << endl;
|
|---|
| 501 |
|
|---|
| 502 | return;
|
|---|
| 503 | }
|
|---|
| 504 |
|
|---|
| 505 | void
|
|---|
| 506 | MonteCarlo::WriteEventData2Csv()
|
|---|
| 507 | {
|
|---|
| 508 | if (mVerbosityLvl > 0) cout << "...writing event data to csv" << endl;
|
|---|
| 509 |
|
|---|
| 510 | for (int i = 0; i < mNumberOfPixels; i++)
|
|---|
| 511 | {
|
|---|
| 512 | WritePixelData2Csv(i);
|
|---|
| 513 | }
|
|---|
| 514 |
|
|---|
| 515 | return;
|
|---|
| 516 | }
|
|---|
| 517 |
|
|---|
| 518 | void
|
|---|
| 519 | MonteCarlo::WriteEventDataNames2Csv()
|
|---|
| 520 | {
|
|---|
| 521 | if (mVerbosityLvl > 0) cout << "...writing event categories to csv" << endl;
|
|---|
| 522 |
|
|---|
| 523 | mCsv << "### [RawData]" << endl;
|
|---|
| 524 | mCsv << "# mEventNumber" << mSeparator;
|
|---|
| 525 | mCsv << "pixelID" << mSeparator;
|
|---|
| 526 | mCsv << "pixelPedestal" << mSeparator;
|
|---|
| 527 | for (int i = 0; i < mNumberOfSamples; i++)
|
|---|
| 528 | {
|
|---|
| 529 | mCsv << "Raw_" << i << mSeparator;
|
|---|
| 530 | }
|
|---|
| 531 | mCsv << endl;
|
|---|
| 532 |
|
|---|
| 533 | return;
|
|---|
| 534 | }
|
|---|
| 535 |
|
|---|
| 536 | void
|
|---|
| 537 | MonteCarlo::WriteEventHeaderNames2Csv()
|
|---|
| 538 | {
|
|---|
| 539 | if (mVerbosityLvl > 0) cout << "...writing event header names to csv" << endl;
|
|---|
| 540 |
|
|---|
| 541 | mCsv << "### [EventHeader]" << endl;
|
|---|
| 542 |
|
|---|
| 543 | mCsv << "# mEventNumber" << mSeparator;
|
|---|
| 544 | mCsv << "mNumberOfBytes" << mSeparator;
|
|---|
| 545 | mCsv << "mIncidentAngle" << mSeparator;
|
|---|
| 546 | mCsv << "mPartId" << mSeparator;
|
|---|
| 547 | mCsv << "mEnergy" << mSeparator;
|
|---|
| 548 | mCsv << "mImpact" << mSeparator;
|
|---|
| 549 | mCsv << "mTelescopePhi" << mSeparator;
|
|---|
| 550 | mCsv << "mTelescopeTheta" << mSeparator;
|
|---|
| 551 | mCsv << "mPhi" << mSeparator;
|
|---|
| 552 | mCsv << "mTheta" << mSeparator;
|
|---|
| 553 | mCsv << "mCorsikaEventNumber" << mSeparator;
|
|---|
| 554 | mCsv << "mPhotElFromShower" << mSeparator;
|
|---|
| 555 | mCsv << "mEvtReuse" << mSeparator;
|
|---|
| 556 | mCsv << "mNumTriggerLvl1" << mSeparator;
|
|---|
| 557 | mCsv << "mFirstInteractionHeight" << mSeparator;
|
|---|
| 558 | mCsv << "mMomentumX" << mSeparator;
|
|---|
| 559 | mCsv << "mMomentumY" << mSeparator;
|
|---|
| 560 | mCsv << "mMomentumZ" << mSeparator;
|
|---|
| 561 | mCsv << "mZd" << mSeparator;
|
|---|
| 562 | mCsv << "mAz" << mSeparator;
|
|---|
| 563 | mCsv << "mX" << mSeparator;
|
|---|
| 564 | mCsv << "mY" << mSeparator;
|
|---|
| 565 | // mCsv << "mWeightedNumPhotons" ;
|
|---|
| 566 | mCsv << endl;
|
|---|
| 567 |
|
|---|
| 568 | return;
|
|---|
| 569 | }
|
|---|
| 570 |
|
|---|
| 571 |
|
|---|
| 572 | void
|
|---|
| 573 | MonteCarlo::WriteEventHeader2Csv()
|
|---|
| 574 | {
|
|---|
| 575 | if (mVerbosityLvl > 0) cout << "...writing event header to csv" << endl;
|
|---|
| 576 |
|
|---|
| 577 | mCsv << mEventNumber << mSeparator;
|
|---|
| 578 | mCsv << mNumberOfBytes << mSeparator;
|
|---|
| 579 | mCsv << mIncidentAngle << mSeparator;
|
|---|
| 580 | mCsv << mPartId << mSeparator;
|
|---|
| 581 | mCsv << mEnergy << mSeparator;
|
|---|
| 582 | mCsv << mImpact << mSeparator;
|
|---|
| 583 | mCsv << mTelescopePhi << mSeparator;
|
|---|
| 584 | mCsv << mTelescopeTheta << mSeparator;
|
|---|
| 585 | mCsv << mPhi << mSeparator;
|
|---|
| 586 | mCsv << mTheta << mSeparator;
|
|---|
| 587 | mCsv << mCorsikaEventNumber << mSeparator;
|
|---|
| 588 | mCsv << mPhotElFromShower << mSeparator;
|
|---|
| 589 | mCsv << mEvtReuse << mSeparator;
|
|---|
| 590 | mCsv << mNumTriggerLvl1 << mSeparator;
|
|---|
| 591 | mCsv << mFirstInteractionHeight << mSeparator;
|
|---|
| 592 | mCsv << mMomentumX << mSeparator;
|
|---|
| 593 | mCsv << mMomentumY << mSeparator;
|
|---|
| 594 | mCsv << mMomentumZ << mSeparator;
|
|---|
| 595 | mCsv << mZd << mSeparator;
|
|---|
| 596 | mCsv << mAz << mSeparator;
|
|---|
| 597 | mCsv << mX << mSeparator;
|
|---|
| 598 | mCsv << mY << mSeparator;
|
|---|
| 599 | // mCsv << mWeightedNumPhotons ;
|
|---|
| 600 | mCsv << endl;
|
|---|
| 601 |
|
|---|
| 602 | return;
|
|---|
| 603 | }
|
|---|
| 604 |
|
|---|
| 605 | void
|
|---|
| 606 | MonteCarlo::WriteRunHeaderNames2Csv()
|
|---|
| 607 | {
|
|---|
| 608 | if (mVerbosityLvl > 0) cout << "...writing run header names to csv" << endl;
|
|---|
| 609 |
|
|---|
| 610 | mCsv << "### [RunHeader]" << endl;
|
|---|
| 611 |
|
|---|
| 612 | mCsv << "# mNumberOfEntries" << mSeparator;
|
|---|
| 613 | mCsv << "mIntendedPulsePos" << mSeparator;
|
|---|
| 614 | // Csv << "mPedestalOffset" << mSeparator;
|
|---|
| 615 | mCsv << "mNumSimulatedShowers" << mSeparator;
|
|---|
| 616 | mCsv << "mNumberOfPixels" << mSeparator;
|
|---|
| 617 | mCsv << "mNumberOfSamples" << mSeparator;
|
|---|
| 618 | // mCsv << "mSampleSize" << mSeparator;
|
|---|
| 619 | mCsv << "mCamDist" << mSeparator;
|
|---|
| 620 | mCsv << "mSourceName" << mSeparator;
|
|---|
| 621 | mCsv << "mSlopeSpectrum" << mSeparator;
|
|---|
| 622 | mCsv << "mEnergyMin" << mSeparator;
|
|---|
| 623 | mCsv << "mEnergyMax" << mSeparator;
|
|---|
| 624 | mCsv << "mZdMin" << mSeparator;
|
|---|
| 625 | mCsv << "mZdMax" << mSeparator;
|
|---|
| 626 | mCsv << "mAzMin" << mSeparator;
|
|---|
| 627 | mCsv << "mAzMax" << mSeparator;
|
|---|
| 628 | mCsv << "mFileNumber" << mSeparator;
|
|---|
| 629 | mCsv << "mRunnumber" << mSeparator;
|
|---|
| 630 | mCsv << "mRunType" << endl;
|
|---|
| 631 |
|
|---|
| 632 | return;
|
|---|
| 633 | }
|
|---|
| 634 |
|
|---|
| 635 |
|
|---|
| 636 | void
|
|---|
| 637 | MonteCarlo::WriteRunHeader2Csv()
|
|---|
| 638 | {
|
|---|
| 639 | if (mVerbosityLvl > 0) cout << "...writing run header to csv" << endl;
|
|---|
| 640 |
|
|---|
| 641 | mCsv << mNumberOfEntries << mSeparator;
|
|---|
| 642 | mCsv << mIntendedPulsePos << mSeparator;
|
|---|
| 643 | // mCsv << mPedestalOffset << mSeparator;
|
|---|
| 644 | mCsv << mNumSimulatedShowers << mSeparator;
|
|---|
| 645 | mCsv << mNumberOfPixels << mSeparator;
|
|---|
| 646 | mCsv << mNumberOfSamples << mSeparator;
|
|---|
| 647 | // mCsv << mSampleSize << mSeparator;
|
|---|
| 648 | mCsv << mCamDist << mSeparator;
|
|---|
| 649 | mCsv << mSourceName << mSeparator;
|
|---|
| 650 | mCsv << mSlopeSpectrum << mSeparator;
|
|---|
| 651 | mCsv << mEnergyMin << mSeparator;
|
|---|
| 652 | mCsv << mEnergyMax << mSeparator;
|
|---|
| 653 | mCsv << mZdMin << mSeparator;
|
|---|
| 654 | mCsv << mZdMax << mSeparator;
|
|---|
| 655 | mCsv << mAzMin << mSeparator;
|
|---|
| 656 | mCsv << mAzMax << mSeparator;
|
|---|
| 657 | mCsv << mFileNumber << mSeparator;
|
|---|
| 658 | mCsv << mRunNumber << mSeparator;
|
|---|
| 659 | mCsv << mRunType << endl;
|
|---|
| 660 |
|
|---|
| 661 | return;
|
|---|
| 662 | }
|
|---|
| 663 |
|
|---|
| 664 | void
|
|---|
| 665 | MonteCarlo::WriteFileHeader2Csv()
|
|---|
| 666 | {
|
|---|
| 667 | if (mVerbosityLvl > 0) cout << "...writing file header to csv" << endl;
|
|---|
| 668 |
|
|---|
| 669 | mCsv << "### ==========================================================="
|
|---|
| 670 | << endl;
|
|---|
| 671 | mCsv << "### = FACT Monte Carlo ="
|
|---|
| 672 | << endl;
|
|---|
| 673 | mCsv << "### ==========================================================="
|
|---|
| 674 | << endl;
|
|---|
| 675 | mCsv << "### = FileInfo: "
|
|---|
| 676 | << endl;
|
|---|
| 677 | mCsv << "### = "
|
|---|
| 678 | << endl;
|
|---|
| 679 | mCsv << "### = Source Name: " << mSourceName << endl;
|
|---|
| 680 | mCsv << "### = Number of Entries: " << mNumberOfEntries << endl;
|
|---|
| 681 | mCsv << "### = Run Number: " << mRunNumber << endl;
|
|---|
| 682 | mCsv << "### = Run Type : " << mRunType << endl;
|
|---|
| 683 | mCsv << "### = File Number: " << mFileNumber << endl;
|
|---|
| 684 | mCsv << "### ==========================================================="
|
|---|
| 685 | << endl ;
|
|---|
| 686 | mCsv << "###"
|
|---|
| 687 | << endl ;
|
|---|
| 688 |
|
|---|
| 689 | return;
|
|---|
| 690 | }
|
|---|
| 691 |
|
|---|
| 692 |
|
|---|
| 693 |
|
|---|
| 694 | //--------------------------------------------------------------------------------
|
|---|
| 695 | //
|
|---|
| 696 | //
|
|---|
| 697 | //
|
|---|
| 698 |
|
|---|
| 699 |
|
|---|
| 700 | // FADC samples (hi gain) of all pixels
|
|---|
| 701 | // This is an array of Byte_t variables. The value of a FADC sample has a
|
|---|
| 702 | // size of n=fNumBytesPerSample bytes. Therefore, a FADC sample value will
|
|---|
| 703 | // occupy n consecutive elements of this array (little endian ordering, i.e,
|
|---|
| 704 | // less significant bits (and bytes) go first.
|
|---|
| 705 | // If m = GetNumHiGainSamples(), the n bytes corresponding to the value of the
|
|---|
| 706 | // i-th FADC sample of the j-th pixel are stored in the n consecutive
|
|---|
| 707 | // positions of this array, starting from fHiGainFadcSamples[j*n*m+i*n]
|
|---|
| 708 |
|
|---|
| 709 |
|
|---|
| 710 | /*
|
|---|
| 711 | int extractMC()
|
|---|
| 712 | {
|
|---|
| 713 |
|
|---|
| 714 | TTree *tree = (TTree*)f.Get("Events");
|
|---|
| 715 |
|
|---|
| 716 | // MArrayB *array = NULL;
|
|---|
| 717 |
|
|---|
| 718 | const float dconv = 2000/4096;
|
|---|
| 719 |
|
|---|
| 720 | unsigned short *short_array;
|
|---|
| 721 | TH1F histo("histo", "histo", 150, 0, 150);
|
|---|
| 722 |
|
|---|
| 723 | TCanvas *canvas = new TCanvas("canvas", "canvas", 0, 0, 400, 400 );
|
|---|
| 724 | canvas->cd();
|
|---|
| 725 | tree->SetBranchAddress("MRawEvtData.", &branch1);
|
|---|
| 726 |
|
|---|
| 727 | tree->GetEntry(22);
|
|---|
| 728 |
|
|---|
| 729 | cout << "Branch1 " << branch1 << endl;
|
|---|
| 730 |
|
|---|
| 731 | array = branch1->GetSamples();
|
|---|
| 732 |
|
|---|
| 733 |
|
|---|
| 734 | for(int i = 0; i < 150; i++){
|
|---|
| 735 | cout << "i_" << i << "\t 1. " << short_array[i] ;
|
|---|
| 736 | short_array[i] = le16toh(short_array[i]);
|
|---|
| 737 |
|
|---|
| 738 | // uint16_t htobe16(uint16_t host_16bits);
|
|---|
| 739 | // uint16_t htole16(uint16_t host_16bits);
|
|---|
| 740 | // uint16_t be16toh(uint16_t big_endian_16bits);
|
|---|
| 741 |
|
|---|
| 742 |
|
|---|
| 743 | cout << "\t 2. " << short_array[i] ;
|
|---|
| 744 | cout << "\t arr. " << (int)array[2*i] << " " << (int)array[2*i+1] << endl;
|
|---|
| 745 |
|
|---|
| 746 | histo.SetBinContent(i, short_array[i]-1000);
|
|---|
| 747 | }
|
|---|
| 748 | histo.DrawCopy();
|
|---|
| 749 | canvas->Modified();
|
|---|
| 750 | canvas->Update();
|
|---|
| 751 |
|
|---|
| 752 | // int m = branch1->GetNumHiGainSamples();
|
|---|
| 753 | // int n = branch1->GetNumBytesPerSample();
|
|---|
| 754 |
|
|---|
| 755 | cout << "GetNumHiGainSamples " << branch1->GetNumHiGainSamples()
|
|---|
| 756 | << " GetNumBytesPerSample " << branch1->GetNumBytesPerSample() << endl;
|
|---|
| 757 |
|
|---|
| 758 | // tree->GetEntry(2);
|
|---|
| 759 |
|
|---|
| 760 | // array = branch1->GetSamples();
|
|---|
| 761 |
|
|---|
| 762 | // cout << "Branch2 " << branch1 << endl;
|
|---|
| 763 |
|
|---|
| 764 | //// for(int i = 1439; i < 1441; i++){
|
|---|
| 765 | //// cout << "Pixel: " << i << endl;
|
|---|
| 766 |
|
|---|
| 767 | // for(int j = 431950; j < 432020; j++){
|
|---|
| 768 | // cout << "Sl: " << j << "\t" ;
|
|---|
| 769 | // cout << (unsigned short) array[j] << "\t";
|
|---|
| 770 | // if ( j%5 == 0 ) cout <<endl;
|
|---|
| 771 |
|
|---|
| 772 | // }
|
|---|
| 773 | // cout << endl;
|
|---|
| 774 | //// }
|
|---|
| 775 |
|
|---|
| 776 | //// tree->SetBranchAddress("MRawEvtData.fHiGainFadcSamples.fArray", array);
|
|---|
| 777 |
|
|---|
| 778 | }
|
|---|
| 779 |
|
|---|
| 780 | */
|
|---|
| 781 |
|
|---|
| 782 |
|
|---|
| 783 |
|
|---|
| 784 |
|
|---|