Index: ct/tools/marsmacros/Mc2Csv =================================================================== --- /fact/tools/marsmacros/Mc2Csv (revision 14605) +++ (revision ) @@ -1,732 +1,0 @@ -#include "MonteCarlo.h" - -MonteCarlo::MonteCarlo() -{ - InitVariables(); - - return; -} - -MonteCarlo::MonteCarlo( - TString filename - ) -{ - InitVariables(); - - mFileName = filename; - - OpenRootFile(); - - LoadEventTree( "Events"); - LoadHeaderTree( "RunHeaders"); - - ReadRunHeader(); - - return; -} - -MonteCarlo::MonteCarlo( - TString filename, - TString evtsTreeName - ) -{ - InitVariables(); - - mFileName = filename; - - OpenRootFile(); - - LoadEventTree( evtsTreeName); - LoadHeaderTree( "RunHeaders"); - - ReadRunHeader(); - - return; -} - -MonteCarlo::MonteCarlo( - TString filename, - TString evtsTreeName, - TString headerTreeName - ) -{ - InitVariables(); - - mFileName = filename; - - OpenRootFile(); - - LoadEventTree( evtsTreeName); - LoadHeaderTree( headerTreeName); - - ReadRunHeader(); - - return; -} - -MonteCarlo::~MonteCarlo(void) -{ - delete[] mpPixel; - - CloseRootFile(); - - delete mpRootFile; - - return; -} - -void -MonteCarlo::InitVariables() -{ - mFileName = ""; - mSeparator = " "; - - // Änderung !!! - mpRootFile = NULL; - // Änderung zuende - - mpEventTree = NULL; - mpHeaderTree = NULL; - mpPixel = NULL; - - // Änderung !!! - mpIntendedPulsePos = NULL; - mpMcRunHeader = NULL; - mpGeomCam = NULL; - mpRawRunHeader = NULL; - mpCorsikaRunHeader = NULL; - - mpElectronicNoise = NULL; - mpRawEventData = NULL; - mpIncidentAngle = NULL; - mpMcEventMetaData = NULL; - mpRawEventHeader = NULL; - mpCorsikaEvtHeader = NULL; - // Änderung zuende - - mpSamples = NULL; - mEventNumber = 0; - mNumberOfEvents = 2; - - mVerbosityLvl = 0; - - return; -} - -void -MonteCarlo::SetVerbosityLevel(int verbLvl) -{ - mVerbosityLvl = verbLvl; - - return; -} - -int -MonteCarlo::GetVerbosityLevel() -{ - return mVerbosityLvl; -} - -void -MonteCarlo::OpenRootFile() -{ - if (mVerbosityLvl > 0) cout << "...opening root file: " << mFileName << endl; - - mpRootFile = new TFile(mFileName, "READ"); - - if (!mpRootFile->IsOpen()) - { - cout << "ERROR - Could not find file " << mFileName << endl; - return ; - } - - return; -} - -void -MonteCarlo::CloseRootFile() -{ - if (mVerbosityLvl > 0) cout << "...closing root file: " << mFileName << endl; - mpRootFile->Close("R"); - - return; -} - -void -MonteCarlo::OpenCsvFile(TString fileName) -{ - mCsvFileName = fileName; - - if (mVerbosityLvl > 0) cout << "...opening csv file: " << mCsvFileName << endl; - - mCsv.open( mCsvFileName ); - - return; -} - -void -MonteCarlo::CloseCsvFile() -{ - if (mVerbosityLvl > 0) cout << "...closing csv file: " << mCsvFileName << endl; - - mCsv.close(); - - return; -} - -void -MonteCarlo::LoadHeaderTree(TString treeName) -{ - if (mVerbosityLvl > 0) cout << "...loading run header tree" << endl; - - mpHeaderTree = (TTree*)mpRootFile->Get(treeName); - - if (mpHeaderTree->IsZombie()) - { - cout << "...could not load tree" << endl; - return; - } - - - - //Set Adresses to Branches in RunHeader-Tree - mpHeaderTree->SetBranchAddress("MGeomCam.", &mpGeomCam); - mpHeaderTree->SetBranchAddress("IntendedPulsePos.", &mpIntendedPulsePos); - mpHeaderTree->SetBranchAddress("MMcRunHeader.", &mpMcRunHeader); - mpHeaderTree->SetBranchAddress("ElectronicNoise.", &mpElectronicNoise); - mpHeaderTree->SetBranchAddress("MRawRunHeader.", &mpRawRunHeader); - mpHeaderTree->SetBranchAddress("MCorsikaRunHeader.",&mpCorsikaRunHeader); - return; -} - -void -MonteCarlo::ReadRunHeader() -{ - if (mVerbosityLvl > 0) - cout << "...reading run header " << mpHeaderTree << endl; - - //Read Values from RunHeader-Tree - - mpHeaderTree->GetEntry(); //This causes problems - - - //Get values from "MGeomCam" Branch - //Getter functions from "MGeomCamFACT.h" - mCamDist = mpGeomCam->GetCameraDist(); - mNumberOfPixels = mpGeomCam->GetNumPixels(); - mNumberOfSectors = mpGeomCam->GetNumSectors(); - mNumberOfAreas = mpGeomCam->GetNumAreas(); - - //Get values from "IntendedPulsePos" Branch - mIntendedPulsePos = mpIntendedPulsePos->GetVal(); - - //Get values from "MMcRunHeader" Branch from event Tree - mNumSimulatedShowers = mpMcRunHeader->GetNumSimulatedShowers(); - - //Get number of Events from event Tree - mNumberOfEntries = mpEventTree->GetEntries(); - if (mVerbosityLvl > 0) - cout << "Event Tree has " << mNumberOfEntries << " entries" << endl; - - //Get values from "MRawRunHeader" Branch - mNumberOfEvents = mpRawRunHeader->GetNumEvents(); // couses problems - mNumberOfEventsRead = mpRawRunHeader->GetNumEventsRead(); - mSamplingFrequency = mpRawRunHeader->GetFreqSampling(); - mSourceName = mpRawRunHeader->GetSourceName(); - mFileNumber = mpRawRunHeader->GetFileNumber(); - mNumberOfSamples = mpRawRunHeader->GetNumSamplesHiGain(); - mNumberOfBytes = mpRawRunHeader->GetNumBytesPerSample(); - mRunNumber = mpRawRunHeader->GetRunNumber(); - mRunType = mpRawRunHeader->GetRunType(); - - - //Get values from "MCorsikaRunHeader" Branch - mSlopeSpectrum = mpCorsikaRunHeader->GetSlopeSpectrum(); - mEnergyMin = mpCorsikaRunHeader->GetEnergyMin(); - mEnergyMax = mpCorsikaRunHeader->GetEnergyMax(); - mZdMin = mpCorsikaRunHeader->GetZdMin(); - mZdMax = mpCorsikaRunHeader->GetZdMax(); - mAzMin = mpCorsikaRunHeader->GetAzMin(); - mAzMax = mpCorsikaRunHeader->GetAzMax(); - - return; -} - -void -MonteCarlo::LoadEventTree(TString treeName) -{ - if (mVerbosityLvl > 0) cout << "...loading event tree" << endl; - - mpEventTree = (TTree*)mpRootFile->Get(treeName); - - if (mpEventTree->IsZombie()) - { - cout << "...could not load tree" << endl; - return; - } - - //Set Adresses to Branches in Events-Tree - mpEventTree->SetBranchAddress("MRawEvtData.", &mpRawEventData); - mpEventTree->SetBranchAddress("IncidentAngle.", &mpIncidentAngle); - mpEventTree->SetBranchAddress("MMcEvt.", &mpMcEventMetaData); - mpEventTree->SetBranchAddress("MRawEvtHeader.", &mpRawEventHeader); - mpEventTree->SetBranchAddress("MCorsikaEvtHeader.", &mpCorsikaEvtHeader); - - return; -} - -void -MonteCarlo::ReadEventMetaData() -{ - if (mVerbosityLvl > 1) cout << "...reading event header" << endl; - - //Get values from "MGeomCamFACT" Branch - mIncidentAngle = mpIncidentAngle->GetVal(); - - //Get values from "MMcEvt" Branch - //The following Getter Functions can be found in MMcEvt.h - mCorsikaEventNumber = mpMcEventMetaData->GetEvtNumber(); - mPhotElFromShower = mpMcEventMetaData->GetPhotElfromShower(); - mPhotElinCamera = mpMcEventMetaData->GetPhotElinCamera(); - //The following Getter Functions can be found in MMcEvtBasic.h - mPartId = mpMcEventMetaData->GetPartId(); - mPartName = mpMcEventMetaData->GetParticleName(mPartId); - mPartSymbol = mpMcEventMetaData->GetParticleSymbol(mPartId); - mEnergy = mpMcEventMetaData->GetEnergy(); - mImpact = mpMcEventMetaData->GetImpact(); - mTelescopePhi = mpMcEventMetaData->GetTelescopePhi(); - mTelescopeTheta = mpMcEventMetaData->GetTelescopeTheta(); - mPhi = mpMcEventMetaData->GetParticlePhi(); - mTheta = mpMcEventMetaData->GetParticleTheta(); - - //Get values from "MRawEventHeader" Branch - mEventNumber = mpRawEventHeader->GetDAQEvtNumber(); - mNumTriggerLvl1 = mpRawEventHeader->GetNumTrigLvl1(); - mNumTriggerLvl2 = mpRawEventHeader->GetNumTrigLvl2(); - - //Get values from "MCorsikaEvtHeader" Branch - mFirstInteractionHeight = mpCorsikaEvtHeader->GetFirstInteractionHeight(); - mEvtReuse = mpCorsikaEvtHeader->GetNumReuse(); - mMomentumX = mpCorsikaEvtHeader->GetMomentum().X(); - mMomentumY = mpCorsikaEvtHeader->GetMomentum().Y(); - mMomentumZ = mpCorsikaEvtHeader->GetMomentum().Z(); - mZd = mpCorsikaEvtHeader->GetZd(); - mAz = mpCorsikaEvtHeader->GetAz(); - mX = mpCorsikaEvtHeader->GetX(); - mY = mpCorsikaEvtHeader->GetY(); - -// mWeightedNumPhotons = mpCorsikaEvtHeader->Get; - // no getter function, no idea how to extract information - - return; -} - -void -MonteCarlo::ReadEventRawData() -{ - if (mVerbosityLvl > 1) cout << "...reading event raw data" << endl; - - // delete Pixel Array before you set refferences for a new one - if (mpPixel != NULL) - { - delete[] mpPixel; - } - mpPixel = new pixel_t[mNumberOfPixels]; - - mpRawEventData->InitRead(mpRawRunHeader); - - int pix_first_sample; -// int pix_last_sample; - - unsigned short* all_raw_data = NULL; - - all_raw_data = (unsigned short*) mpRawEventData->GetSamples(); - - if (mVerbosityLvl > 1) - cout << "...pixel progress: "; - - for ( int i = 0; i < mNumberOfPixels; i++ ) - { - if (mVerbosityLvl > 1){ - if ( !(i%(mNumberOfPixels/10) )) - cout << i/(mNumberOfPixels*1.0)*100 << "%.." ; - if ( i == mNumberOfPixels-1) - cout << "100% " ; - } - - - mpPixel[i].pedestal = mpElectronicNoise[0][i].GetPedestal(); - mpPixel[i].SoftId = mpRawEventData->GetPixelIds()[i]; - - pix_first_sample = i*mNumberOfSamples; - - // point beginning of pixel's rawdata array to address - // of pixel's first sample's adress - - mpPixel[i].rawData = &(all_raw_data[pix_first_sample]); - -// for (int j = pix_first_sample; j < pix_last_sample; j++) -// { -// cout << mpPixel[i].rawData[j] << " "; -// } -// cout << endl; - -// // this could be improved by -// DONE BUT NOT TESTED -// for (int j = i*mSampleSize; j < (i+1)*mSampleSize; j++) -// { -// } - } - - if (mVerbosityLvl > 1) - cout << endl; - - return; -} - -void -MonteCarlo::ReadEvent(int Event) -{ - if (mVerbosityLvl > 0) cout << endl - << "====================" << endl - << "...reading Event: " << Event << endl - << "====================" << endl; - - - //load certain event from tree - mpEventTree->GetEntry(Event); - - ReadEventMetaData(); - ReadEventRawData(); - - return; -} - - -void -MonteCarlo::WriteMc2Csv(TString filename) -{ - cout << "...writing mc to csv: " << filename << endl; - cout << "...processing " << mNumberOfEntries << "Events" << endl; - - mCsvFileName = filename; - OpenCsvFile(mCsvFileName); - - WriteFileHeader2Csv(); - WriteRunHeaderNames2Csv(); - WriteRunHeader2Csv(); - WriteEventHeaderNames2Csv(); - WriteEventDataNames2Csv(); - -// for (int evt = 0; evt < mNumberOfEvents; evt++) - for (int evt = 0; evt < mNumberOfEntries; evt++) - { - ReadEvent(evt); - WriteEventHeader2Csv(); - WriteEventData2Csv(); - } - - cout << endl << "...conversion done " << endl; - - CloseCsvFile(); - return; -} - - -void -MonteCarlo::WritePixelData2Csv(int pixelID) -{ - if (mVerbosityLvl > 3) cout << "...writing pixel data to csv" << endl; - mCsv << mEventNumber << mSeparator; - mCsv << mpPixel[pixelID].SoftId << mSeparator; - mCsv << mpPixel[pixelID].pedestal << mSeparator; - - for (int i = 0; i < mNumberOfSamples; i++) - { - mCsv << mpPixel[pixelID].rawData[i] << mSeparator; - } - mCsv << endl; - - return; -} - -void -MonteCarlo::WriteEventData2Csv() -{ - if (mVerbosityLvl > 0) cout << "...writing event data to csv" << endl; - - for (int i = 0; i < mNumberOfPixels; i++) - { - WritePixelData2Csv(i); - } - - return; -} - -void -MonteCarlo::WriteEventDataNames2Csv() -{ - if (mVerbosityLvl > 0) cout << "...writing event categories to csv" << endl; - - mCsv << "### [RawData]" << endl; - mCsv << "# mEventNumber" << mSeparator; - mCsv << "pixelID" << mSeparator; - mCsv << "pixelPedestal" << mSeparator; - for (int i = 0; i < mNumberOfSamples; i++) - { - mCsv << "Raw_" << i << mSeparator; - } - mCsv << endl; - - return; -} - -void -MonteCarlo::WriteEventHeaderNames2Csv() -{ - if (mVerbosityLvl > 0) cout << "...writing event header names to csv" << endl; - - mCsv << "### [EventHeader]" << endl; - - mCsv << "# mEventNumber" << mSeparator; - mCsv << "mNumberOfBytes" << mSeparator; - mCsv << "mIncidentAngle" << mSeparator; - mCsv << "mPartId" << mSeparator; - mCsv << "mEnergy" << mSeparator; - mCsv << "mImpact" << mSeparator; - mCsv << "mTelescopePhi" << mSeparator; - mCsv << "mTelescopeTheta" << mSeparator; - mCsv << "mPhi" << mSeparator; - mCsv << "mTheta" << mSeparator; - mCsv << "mCorsikaEventNumber" << mSeparator; - mCsv << "mPhotElFromShower" << mSeparator; - mCsv << "mEvtReuse" << mSeparator; - mCsv << "mNumTriggerLvl1" << mSeparator; - mCsv << "mFirstInteractionHeight" << mSeparator; - mCsv << "mMomentumX" << mSeparator; - mCsv << "mMomentumY" << mSeparator; - mCsv << "mMomentumZ" << mSeparator; - mCsv << "mZd" << mSeparator; - mCsv << "mAz" << mSeparator; - mCsv << "mX" << mSeparator; - mCsv << "mY" << mSeparator; -// mCsv << "mWeightedNumPhotons" ; - mCsv << endl; - - return; -} - - -void -MonteCarlo::WriteEventHeader2Csv() -{ - if (mVerbosityLvl > 0) cout << "...writing event header to csv" << endl; - - mCsv << mEventNumber << mSeparator; - mCsv << mNumberOfBytes << mSeparator; - mCsv << mIncidentAngle << mSeparator; - mCsv << mPartId << mSeparator; - mCsv << mEnergy << mSeparator; - mCsv << mImpact << mSeparator; - mCsv << mTelescopePhi << mSeparator; - mCsv << mTelescopeTheta << mSeparator; - mCsv << mPhi << mSeparator; - mCsv << mTheta << mSeparator; - mCsv << mCorsikaEventNumber << mSeparator; - mCsv << mPhotElFromShower << mSeparator; - mCsv << mEvtReuse << mSeparator; - mCsv << mNumTriggerLvl1 << mSeparator; - mCsv << mFirstInteractionHeight << mSeparator; - mCsv << mMomentumX << mSeparator; - mCsv << mMomentumY << mSeparator; - mCsv << mMomentumZ << mSeparator; - mCsv << mZd << mSeparator; - mCsv << mAz << mSeparator; - mCsv << mX << mSeparator; - mCsv << mY << mSeparator; -// mCsv << mWeightedNumPhotons ; - mCsv << endl; - - return; -} - -void -MonteCarlo::WriteRunHeaderNames2Csv() -{ - if (mVerbosityLvl > 0) cout << "...writing run header names to csv" << endl; - - mCsv << "### [RunHeader]" << endl; - - mCsv << "# mNumberOfEntries" << mSeparator; - mCsv << "mIntendedPulsePos" << mSeparator; -// Csv << "mPedestalOffset" << mSeparator; - mCsv << "mNumSimulatedShowers" << mSeparator; - mCsv << "mNumberOfPixels" << mSeparator; - mCsv << "mNumberOfSamples" << mSeparator; -// mCsv << "mSampleSize" << mSeparator; - mCsv << "mCamDist" << mSeparator; - mCsv << "mSourceName" << mSeparator; - mCsv << "mSlopeSpectrum" << mSeparator; - mCsv << "mEnergyMin" << mSeparator; - mCsv << "mEnergyMax" << mSeparator; - mCsv << "mZdMin" << mSeparator; - mCsv << "mZdMax" << mSeparator; - mCsv << "mAzMin" << mSeparator; - mCsv << "mAzMax" << mSeparator; - mCsv << "mFileNumber" << mSeparator; - mCsv << "mRunnumber" << mSeparator; - mCsv << "mRunType" << endl; - - return; -} - - -void -MonteCarlo::WriteRunHeader2Csv() -{ - if (mVerbosityLvl > 0) cout << "...writing run header to csv" << endl; - - mCsv << mNumberOfEntries << mSeparator; - mCsv << mIntendedPulsePos << mSeparator; -// mCsv << mPedestalOffset << mSeparator; - mCsv << mNumSimulatedShowers << mSeparator; - mCsv << mNumberOfPixels << mSeparator; - mCsv << mNumberOfSamples << mSeparator; -// mCsv << mSampleSize << mSeparator; - mCsv << mCamDist << mSeparator; - mCsv << mSourceName << mSeparator; - mCsv << mSlopeSpectrum << mSeparator; - mCsv << mEnergyMin << mSeparator; - mCsv << mEnergyMax << mSeparator; - mCsv << mZdMin << mSeparator; - mCsv << mZdMax << mSeparator; - mCsv << mAzMin << mSeparator; - mCsv << mAzMax << mSeparator; - mCsv << mFileNumber << mSeparator; - mCsv << mRunNumber << mSeparator; - mCsv << mRunType << endl; - - return; -} - -void -MonteCarlo::WriteFileHeader2Csv() -{ - if (mVerbosityLvl > 0) cout << "...writing file header to csv" << endl; - - mCsv << "### ===========================================================" - << endl; - mCsv << "### = FACT Monte Carlo =" - << endl; - mCsv << "### ===========================================================" - << endl; - mCsv << "### = FileInfo: " - << endl; - mCsv << "### = " - << endl; - mCsv << "### = Source Name: " << mSourceName << endl; - mCsv << "### = Number of Entries: " << mNumberOfEntries << endl; - mCsv << "### = Run Number: " << mRunNumber << endl; - mCsv << "### = Run Type : " << mRunType << endl; - mCsv << "### = File Number: " << mFileNumber << endl; - mCsv << "### ===========================================================" - << endl ; - mCsv << "###" - << endl ; - - return; -} - - - -//-------------------------------------------------------------------------------- -// -// -// - - -// FADC samples (hi gain) of all pixels -// This is an array of Byte_t variables. The value of a FADC sample has a -// size of n=fNumBytesPerSample bytes. Therefore, a FADC sample value will -// occupy n consecutive elements of this array (little endian ordering, i.e, -// less significant bits (and bytes) go first. -// If m = GetNumHiGainSamples(), the n bytes corresponding to the value of the -// i-th FADC sample of the j-th pixel are stored in the n consecutive -// positions of this array, starting from fHiGainFadcSamples[j*n*m+i*n] - - -/* -int extractMC() -{ - - TTree *tree = (TTree*)f.Get("Events"); - -// MArrayB *array = NULL; - - const float dconv = 2000/4096; - - unsigned short *short_array; - TH1F histo("histo", "histo", 150, 0, 150); - - TCanvas *canvas = new TCanvas("canvas", "canvas", 0, 0, 400, 400 ); - canvas->cd(); - tree->SetBranchAddress("MRawEvtData.", &branch1); - - tree->GetEntry(22); - - cout << "Branch1 " << branch1 << endl; - - array = branch1->GetSamples(); - - - for(int i = 0; i < 150; i++){ - cout << "i_" << i << "\t 1. " << short_array[i] ; - short_array[i] = le16toh(short_array[i]); - -// uint16_t htobe16(uint16_t host_16bits); -// uint16_t htole16(uint16_t host_16bits); -// uint16_t be16toh(uint16_t big_endian_16bits); - - - cout << "\t 2. " << short_array[i] ; - cout << "\t arr. " << (int)array[2*i] << " " << (int)array[2*i+1] << endl; - - histo.SetBinContent(i, short_array[i]-1000); - } - histo.DrawCopy(); - canvas->Modified(); - canvas->Update(); - -// int m = branch1->GetNumHiGainSamples(); -// int n = branch1->GetNumBytesPerSample(); - - cout << "GetNumHiGainSamples " << branch1->GetNumHiGainSamples() - << " GetNumBytesPerSample " << branch1->GetNumBytesPerSample() << endl; - -// tree->GetEntry(2); - -// array = branch1->GetSamples(); - -// cout << "Branch2 " << branch1 << endl; - -//// for(int i = 1439; i < 1441; i++){ -//// cout << "Pixel: " << i << endl; - -// for(int j = 431950; j < 432020; j++){ -// cout << "Sl: " << j << "\t" ; -// cout << (unsigned short) array[j] << "\t"; -// if ( j%5 == 0 ) cout <SetBranchAddress("MRawEvtData.fHiGainFadcSamples.fArray", array); - -} - -*/ - - - -