diff --git a/Makefile b/Makefile index 4c72ba360..2e9c90b1e 100644 --- a/Makefile +++ b/Makefile @@ -286,6 +286,7 @@ EVNOBJECTS = ./obj/VVirtualDataReader.o \ ./obj/VSpectralWeight.o ./obj/VSpectralWeight_Dict.o \ ./obj/VTableLookupRunParameter.o ./obj/VTableLookupRunParameter_Dict.o \ ./obj/VPedestalCalculator.o \ + ./obj/VIPRCalculator.o \ ./obj/VDeadChannelFinder.o \ ./obj/VSpecialChannel.o \ ./obj/VDeadTime.o \ diff --git a/inc/VCalibrationData.h b/inc/VCalibrationData.h index e7917d3a2..718b66b73 100644 --- a/inc/VCalibrationData.h +++ b/inc/VCalibrationData.h @@ -143,7 +143,7 @@ class VCalibrationData // IPR graphs map< unsigned int, TGraphErrors* > fGraphIPRGraph; // one IPR graph per summation window - + map< unsigned int, TGraphErrors* > fGraphTSIPRGraph; // one IPR graph per time slice VCalibrationData( unsigned int iTel, string iDir, string iPedfile, string iGainfile, string iTofffile, string iPedLowGainfile, string iGainLowGainFile = "", string iToffLowGainFile = "", string iLowGainMultFile = "", @@ -237,6 +237,7 @@ class VCalibrationData } } TGraphErrors* getIPRGraph( unsigned int iSumWindow, bool iMakeNewGraph = false ); + TGraphErrors* getIPRGraphTimeSlice( bool iMakeNewGraph = false, unsigned int TimeSlice = 0 ); TH1F* getLowGainMultiplierDistribution() { return getHistoDist( C_LOWGAIN, true ); diff --git a/inc/VCalibrator.h b/inc/VCalibrator.h index 0ed738ef2..4f1effaa9 100644 --- a/inc/VCalibrator.h +++ b/inc/VCalibrator.h @@ -5,6 +5,7 @@ #include "VImageBaseAnalyzer.h" #include "VPedestalCalculator.h" +#include "VIPRCalculator.h" #include "VDB_CalibrationInfo.h" #include "VSQLTextFileReader.h" @@ -121,7 +122,7 @@ class VCalibrator : public VImageBaseAnalyzer void setCalibrationFileNames(); void writeGains( bool iLowGain = false ); - void writePeds( bool iLowGain, VPedestalCalculator* iP = 0, bool iWriteAsciiFile = true ); + void writePeds( bool iLowGain, VPedestalCalculator* iP = 0, bool iWriteAsciiFile = true, VIPRCalculator* fIPRCalculator = 0 ); void writeTOffsets( bool iLowGain = false ); void writeAverageTZeros( bool iLowGain = false ); bool writeIPRgraphs( string iFile = "" ); @@ -135,7 +136,7 @@ class VCalibrator : public VImageBaseAnalyzer void calculatePedestals( bool iLowGain = false ); void calculateGainsAndTOffsets( bool iLowGain = false ); unsigned int getNumberOfEventsUsedInCalibration( int iTelID, int iType ); - void initialize(); - void terminate( VPedestalCalculator* ); + void initialize(VIPRCalculator *fIPRCalculator ); + void terminate( VPedestalCalculator*, VIPRCalculator* ); }; #endif diff --git a/inc/VEventLoop.h b/inc/VEventLoop.h index e8562e014..168cbdcdb 100644 --- a/inc/VEventLoop.h +++ b/inc/VEventLoop.h @@ -21,6 +21,7 @@ #include "VBFDataReader.h" #endif #include "VPedestalCalculator.h" +#include "VIPRCalculator.h" #include "VEvndispRunParameter.h" #include "VDeadPixelOrganizer.h" @@ -40,6 +41,7 @@ class VEventLoop : public VEvndispData enum E_runmode {R_ANA, R_PED, R_GTO, R_BCK, R_DST, R_GTOLOW, R_PEDLOW, R_TZERO, R_TZEROLOW }; VCalibrator* fCalibrator; //!< default calibration class + VIPRCalculator* fIPRCalculator; VPedestalCalculator* fPedestalCalculator; //!< default pedestal calculator VImageAnalyzer* fAnalyzer; //!< default analyzer class VArrayAnalyzer* fArrayAnalyzer; //!< default array analyzer @@ -55,6 +57,7 @@ class VEventLoop : public VEvndispData bool fAnalyzeMode; //!< used for gotoEvent (go through file without analyse events) bool bMCSetAtmosphericID; vector< bool > fBoolPrintSample; + bool fIPRTimeSlices; bool fCutTelescope; //!< cuts apply only to one telescope int fNCutNArrayTrigger; //!< show only events with more than fNCutNArrayTrigger triggered telescopes diff --git a/inc/VEvndispData.h b/inc/VEvndispData.h index be1574d03..954e23942 100644 --- a/inc/VEvndispData.h +++ b/inc/VEvndispData.h @@ -523,6 +523,14 @@ class VEvndispData { return fCalData[fTelID]->getIPRGraph( iSumWindow, iMakeNewGraph ); } + TGraphErrors* getIPRGraphTimeSlice( unsigned int TimeSlice = 0 ) + { + return fCalData[fTelID]->getIPRGraphTimeSlice( false, TimeSlice ); + } + TGraphErrors* getIPRGraphTimeSlice( bool iMakeNewGraph = false, unsigned int TimeSlice = 0 ) + { + return fCalData[fTelID]->getIPRGraphTimeSlice( iMakeNewGraph, TimeSlice ); + } float getL1Rate( unsigned int iChannel ) { if( fDB_PixelDataReader ) diff --git a/inc/VEvndispRunParameter.h b/inc/VEvndispRunParameter.h index fb9a4aaa8..07d28bbd8 100644 --- a/inc/VEvndispRunParameter.h +++ b/inc/VEvndispRunParameter.h @@ -225,6 +225,7 @@ class VEvndispRunParameter : public TNamed, public VGlobalRunParameter bool ifWriteGraphsToFile; // flag to write run-time NN image cleaning settings to root file (read from NN image cleaning input card) bool ifReadIPRfromDatabase; // flag to read IPR from external IPR database bool ifCreateIPRdatabase; // flag to create external IPR database + bool ifReadIPRfromDSTFile; // flat to read IPR graph from DST file string freconstructionparameterfile; // reconstruction parameter file // MC parameters diff --git a/inc/VIPRCalculator.h b/inc/VIPRCalculator.h new file mode 100644 index 000000000..39e8d3e16 --- /dev/null +++ b/inc/VIPRCalculator.h @@ -0,0 +1,57 @@ +#ifndef VIPRCALCULATOR_H +#define VIPRCALCULATOR_H + +#include +#include +#include + +#include "TClonesArray.h" +#include "TFile.h" +#include "TH1F.h" +#include "TLeaf.h" +#include "TMath.h" +#include "TProfile.h" +#include "TSystem.h" +#include "TTree.h" + +#include +#include +#include +#include +#include + +using namespace std; + +class VIPRCalculator : public VImageBaseAnalyzer +{ + + private: + VEvndispData* fData; + vector< string > fPedFileNameC; + bool fIPRTimeSlices; + bool fIPRAverageTel; // flag to make average of all telescopes IPR in case there is not enough statistics to produce IPR graphs + bool fIPRInTimeSlices; + int fPedPerTelescopeTypeMinCnt; + TH1F* FillIPRHistogram( unsigned int iSummationWindow, unsigned int i_tel); + void definePedestalFile( std::vector fPedFileNameCalibrator ); + TH1F* initializeIPRHistogram( unsigned int iSummationWindow, unsigned int i_tel); + bool copyIPRTelAveraged( unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel ); + TH1F* calculateIPRGraphAveraged( unsigned int iSummationWindow ); + public: + vector>>> fpedcal_histo_storage; + bool calculateIPRGraphs( std::vector fPedFileNameCalibrator ); + bool calculateIPRGraphs( string iPedFileName, unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel ); + bool calculateIPRGraphsTimeSlices( string iPedFileName, int TS, unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel ); + bool writeIPRgraphs( map>> &hped_vec, string iFile = "" ); + void fillIPRPedestalHisto(const int telID, const int NTimeSlices,const vector>>& fpedcal_histo ); + void fillIPRPedestalHisto(); + TH1F* getIPRPedestalHisto(const int telID, const int ts, const int pixel, const int sw); + bool clearHistos(); + vector>>> getStorageHist(); + + VIPRCalculator(); + ~VIPRCalculator() {} + + void initialize(); +}; +#endif diff --git a/inc/VPedestalCalculator.h b/inc/VPedestalCalculator.h index 5bf71b24b..a4c077578 100644 --- a/inc/VPedestalCalculator.h +++ b/inc/VPedestalCalculator.h @@ -6,6 +6,7 @@ #include "VImageBaseAnalyzer.h" #include "VGlobalRunParameter.h" #include "VSkyCoordinatesUtilities.h" +#include "VIPRCalculator.h" #include "TDirectory.h" #include "TFile.h" @@ -47,6 +48,7 @@ class VPedestalCalculator : public VImageBaseAnalyzer vector< vector< vector< float > > > fpedcal_mean; vector< vector< vector< float > > > fpedcal_mean2; vector< vector< vector< TH1F* > > > fpedcal_histo; + vector< vector< vector< TH1F* > > > fpedcal_histo_slidingw; vector< vector< float > > v_temp_pedEntries; vector< vector< float > > v_temp_ped; @@ -62,7 +64,7 @@ class VPedestalCalculator : public VImageBaseAnalyzer void reset(); public: - + vector< int > NTimeSlices; vector< vector< int > > v_MJD; //! [telid][time slice] vector< vector< double > > v_time; //! [telid][time slice] //! [telid][time slice][npixel][summation window] @@ -77,7 +79,7 @@ class VPedestalCalculator : public VImageBaseAnalyzer VPedestalCalculator(); ~VPedestalCalculator() {} - void doAnalysis( bool iLowGain = false ); + void doAnalysis( bool iLowGain = false, VIPRCalculator *fIPRCalculator = 0); vector< TTree* > getPedestalTree() { return fTree; @@ -85,6 +87,6 @@ class VPedestalCalculator : public VImageBaseAnalyzer bool initialize(); bool initialize( bool ibCalibrationRun, unsigned int iNPixel, double iLengthofTimeSlice, int iSumFirst, int iSumWindow, double iRunStartTime = -99., double iRunStoppTime = -99. ); - void terminate( bool iWrite = true, bool bDebug_IO = false ); + void terminate( bool iWrite = true, bool bDebug_IO = false, VIPRCalculator* fIPRCalculator = 0); }; #endif diff --git a/src/VCalibrationData.cpp b/src/VCalibrationData.cpp index 5f7d7ee5c..1badb4c2c 100644 --- a/src/VCalibrationData.cpp +++ b/src/VCalibrationData.cpp @@ -1043,6 +1043,25 @@ TGraphErrors* VCalibrationData::getIPRGraph( unsigned int iSumWindow, bool iMake return 0; } +TGraphErrors* VCalibrationData::getIPRGraphTimeSlice( bool iMakeNewGraph, unsigned int TimeSlice ) +{ + + if( fGraphTSIPRGraph.find( TimeSlice ) != fGraphTSIPRGraph.end() && fGraphTSIPRGraph[TimeSlice] ) + { + return fGraphTSIPRGraph[TimeSlice]; + } + else if( iMakeNewGraph ) + { + fGraphTSIPRGraph[TimeSlice] = new TGraphErrors( 1 ); + fGraphTSIPRGraph[TimeSlice]->SetTitle( "" ); + char hname[200]; + sprintf( hname, "IRPFGraph_TelID%d_TimeSlice%d", fTelID, TimeSlice ); + fGraphTSIPRGraph[TimeSlice]->SetName( hname ); + return fGraphTSIPRGraph[TimeSlice]; + } + return 0; +} + void VCalibrationData::setIPRGraph( unsigned int iSumWindow, TGraphErrors* g ) { fGraphIPRGraph[iSumWindow] = g; diff --git a/src/VCalibrator.cpp b/src/VCalibrator.cpp index fcda4786a..7b9a88e1f 100644 --- a/src/VCalibrator.cpp +++ b/src/VCalibrator.cpp @@ -359,7 +359,8 @@ void VCalibrator::calculatePedestals( bool iLowGain ) this might be an ascii file and/or a root file */ -void VCalibrator::writePeds( bool iLowGain, VPedestalCalculator* iPedestalCalculator, bool iWriteAsciiFile ) +//MK test +void VCalibrator::writePeds( bool iLowGain, VPedestalCalculator* iPedestalCalculator, bool iWriteAsciiFile, VIPRCalculator* fIPRCalculator ) { if( getDebugFlag() ) { @@ -484,6 +485,7 @@ void VCalibrator::writePeds( bool iLowGain, VPedestalCalculator* iPedestalCalcul std::ostringstream iSname; iSname << "distributions_" << telType; TDirectory* i_dist = getPedestalRootFile( telType )->mkdir( iSname.str().c_str() ); + if( i_dist->cd() ) { i_dist->cd(); @@ -497,6 +499,7 @@ void VCalibrator::writePeds( bool iLowGain, VPedestalCalculator* iPedestalCalcul } } } + for( unsigned int i = 0; i < hpedPerTelescopeType[telType].size(); i++ ) { for( unsigned int j = 0; j < hpedPerTelescopeType[telType][i].size(); j++ ) @@ -507,11 +510,27 @@ void VCalibrator::writePeds( bool iLowGain, VPedestalCalculator* iPedestalCalcul } } } + //cout << "MK sum window: " << getSumWindow() << " " << getRunParameter()->fCalibrationSumWindow << endl; + for (unsigned int ts = 0; ts < (fIPRCalculator->getStorageHist())[tel].size(); ts++) + { + for (unsigned int p = 0; p < getNChannels(); p++) + { + int sw = getSumWindow() - 6; + if ( fIPRCalculator->getIPRPedestalHisto(tel, ts, p, sw) ) + { + fIPRCalculator->getIPRPedestalHisto(tel, ts, p, sw)->SetName(Form( "hpedTimeSlices_Tel%d_TS%d_Pix%d_SW%d", (int)telType, ts, p, sw + 6 ) ); + fIPRCalculator->getIPRPedestalHisto(tel, ts, p, sw)->SetTitle( Form("hpedTimeSlices_Tel%d_TS%d_Pix%d_SW%d", (int)telType, ts, p, sw + 6 ) ); + fIPRCalculator->getIPRPedestalHisto(tel, ts, p, sw)->Write(); + } + } + } } iFileWritten[telType] = true; } } // end loop over all telescopes // delete all histograms + cout << "MK clearing histos" << endl; + fIPRCalculator->clearHistos(); map< ULong64_t, TClonesArray* >::iterator i_PedestalsHistoClonesArray_iter; for( i_PedestalsHistoClonesArray_iter = fPedestalsHistoClonesArray.begin(); i_PedestalsHistoClonesArray_iter != fPedestalsHistoClonesArray.end(); i_PedestalsHistoClonesArray_iter++ ) @@ -584,7 +603,7 @@ void VCalibrator::writePeds( bool iLowGain, VPedestalCalculator* iPedestalCalcul } } } - writeIPRgraphs( fPedSingleOutFile->GetName() ); + fIPRCalculator->writeIPRgraphs( hped_vec, fPedSingleOutFile->GetName() ); } } @@ -1589,11 +1608,12 @@ void VCalibrator::writeTOffsets( bool iLowGain ) } -void VCalibrator::terminate( VPedestalCalculator* iP ) +void VCalibrator::terminate( VPedestalCalculator* iP, VIPRCalculator* fIPRCalculator ) { if( fRunPar->frunmode == 1 || fRunPar->frunmode == 6 ) { - writePeds( fRunPar->frunmode == 6, iP, !fRunPar->fPedestalSingleRootFile ); + writePeds( fRunPar->frunmode == 6, iP, !fRunPar->fPedestalSingleRootFile, fIPRCalculator ); + } else if( fRunPar->frunmode == 2 || fRunPar->frunmode == 5 ) { @@ -3090,8 +3110,7 @@ void VCalibrator::readTOffsets( bool iLowGain ) } - -void VCalibrator::initialize() +void VCalibrator::initialize( VIPRCalculator* fIPRCalculator ) { if( fDebug ) { @@ -3221,7 +3240,7 @@ void VCalibrator::initialize() // if needed: write IPR graphs to disk if( getRunParameter()->ifCreateIPRdatabase == true && getRunParameter()->ifReadIPRfromDatabase == false ) { - writeIPRgraphs(); + fIPRCalculator->writeIPRgraphs(hped_vec, ""); } // initialize dead channel finder @@ -3237,7 +3256,7 @@ void VCalibrator::initialize() && fRunPar->frunmode != 1 && fRunPar->frunmode != 6 ) { - calculateIPRGraphs(); + fIPRCalculator->calculateIPRGraphs( fPedFileNameC ); } } } @@ -4376,7 +4395,16 @@ bool VCalibrator::readCalibrationDatafromDSTFiles( string iDSTfile ) //////////////////////////////////////////////////////////////////////////// // read IPR graph from dst root file (for direct usage or creation of database ) - if( getRunParameter()->ifReadIPRfromDatabase == true || getRunParameter()->ifCreateIPRdatabase == true ) + if( getRunParameter()->ifReadIPRfromDSTFile == true ) + { + cout << "\t reading IPR graphs for NN image cleaning from DST file" << endl; + for( int i = 0; i < t->GetEntries(); i++ ) + { + setTelID( i ); + readIPRGraph_from_DSTFile( iDSTfile, getSumWindow(), getTelType( i ) ); + } + } + else if( getRunParameter()->ifReadIPRfromDatabase == true || getRunParameter()->ifCreateIPRdatabase == true ) { cout << "\t reading IPR graphs for NN image cleaning" << endl; for( int i = 0; i < t->GetEntries(); i++ ) diff --git a/src/VEventLoop.cpp b/src/VEventLoop.cpp index 0931cb181..5e8b320d6 100644 --- a/src/VEventLoop.cpp +++ b/src/VEventLoop.cpp @@ -101,6 +101,7 @@ VEventLoop::VEventLoop( VEvndispRunParameter* irunparameter ) fCalibrated.push_back( false ); } fCalibrator = new VCalibrator(); + fIPRCalculator = new VIPRCalculator(); // create data summarizer fDST = new VDST( ( fRunMode == R_DST ), ( fRunPar->fsourcetype == 1 || fRunPar->fsourcetype == 2 || fRunPar->fsourcetype == 6 ) ); @@ -475,11 +476,15 @@ bool VEventLoop::initEventLoop( string iFileName ) // initialize analyzers (output files are created as well here) initializeAnalyzers(); - + if( fIPRCalculator ) + { + cout << "initializing IPR calculator" << endl; + fIPRCalculator->initialize(); + } // create calibrators, analyzers, etc. at first event if( fCalibrator ) { - fCalibrator->initialize(); + fCalibrator->initialize( fIPRCalculator ); } // initialize pedestal calculator @@ -756,7 +761,7 @@ void VEventLoop::shutdown() // write pedestal variation calculations to output file if( ( fRunPar->frunmode == R_ANA ) && fRunPar->fPedestalsInTimeSlices && fPedestalCalculator ) { - fPedestalCalculator->terminate( true, fDebug_writing ); + fPedestalCalculator->terminate( true, fDebug_writing, fIPRCalculator ); } // calculate and write deadtime calculation to disk if( fDeadTime && !isMC() ) @@ -817,11 +822,11 @@ void VEventLoop::shutdown() if( fRunPar->frunmode == R_PED && fRunPar->fPedestalsInTimeSlices && fPedestalCalculator ) { iP = fPedestalCalculator; - fPedestalCalculator->terminate( false ); + fPedestalCalculator->terminate( false, false, fIPRCalculator ); } if( fCalibrator ) { - fCalibrator->terminate( iP ); + fCalibrator->terminate( iP , fIPRCalculator ); } } // write data summary @@ -1292,6 +1297,7 @@ int VEventLoop::analyzeEvent() // analysis case R_ANA: // analysis mode // ignore pedestal events (important for VBF only) + cout << getEventNumber() << endl; #ifndef NOVBF if( fReader->getATEventType() != VEventType::PED_TRIGGER ) #endif @@ -1452,7 +1458,7 @@ int VEventLoop::analyzeEvent() setTelID( fRunPar->fTelToAnalyze[i] ); if( fRunPar->fPedestalsInTimeSlices && !fReader->wasLossyCompressed() ) { - fPedestalCalculator->doAnalysis( fRunMode == R_PEDLOW ); + fPedestalCalculator->doAnalysis( fRunMode == R_PEDLOW, fIPRCalculator ); } } } diff --git a/src/VEvndispRunParameter.cpp b/src/VEvndispRunParameter.cpp index 75b4746d1..fbb148653 100644 --- a/src/VEvndispRunParameter.cpp +++ b/src/VEvndispRunParameter.cpp @@ -270,6 +270,7 @@ VEvndispRunParameter::VEvndispRunParameter( bool bSetGlobalParameter ) : VGlobal ifWriteGraphsToFile = false; ifReadIPRfromDatabase = false; ifCreateIPRdatabase = false; + ifReadIPRfromDSTFile = false; fNNGraphsFile = "IPRgraph.root"; fIPRdatabase = ""; diff --git a/src/VIPRCalculator.cpp b/src/VIPRCalculator.cpp new file mode 100644 index 000000000..3ce56d6fa --- /dev/null +++ b/src/VIPRCalculator.cpp @@ -0,0 +1,839 @@ +/*! \class IPRCalculator + * \brief calculation and save IPR GRAPH + * + * */ + +#include "VIPRCalculator.h" + +VIPRCalculator::VIPRCalculator() +{ + fIPRTimeSlices = true; + fIPRAverageTel = false; + fPedPerTelescopeTypeMinCnt = 1.e3; //1.E5; // minimal counter for IPR measurements + fpedcal_histo_storage.resize( getTeltoAna().size() ); +} + +void VIPRCalculator::initialize() +{ + if( fDebug ) + { + cout << "VIPRCalculator::initialize()" << endl; + } + + // set the data readers + initializeDataReader(); + if( !initializeDataReader() ) + { + cout << "VIPRCalculator::initialize, error: cannot initialize data readers" << endl; + cout << "exiting..." << endl; + exit( EXIT_FAILURE ); + } + +} + +void VIPRCalculator::definePedestalFile( std::vector fPedFileNameCalibrator ) +{ + for( unsigned int i = 0; i < fPedFileNameCalibrator.size() ; i++ ) + { + fPedFileNameC.push_back( fPedFileNameCalibrator[i] ); + } +} + +bool VIPRCalculator::calculateIPRGraphs( std::vector fPedFileNameCalibrator ) +{ + + definePedestalFile( fPedFileNameCalibrator ); + for( unsigned int i = 0; i < getTeltoAna().size(); i++ ) + { + setTelID( getTeltoAna()[i] ); + + // first find dead channels + // (ignore low gains here) + findDeadChans( false ); + // calculate IPR graphs + if( fIPRAverageTel == false ) + { + calculateIPRGraphs( fPedFileNameC[getTeltoAna()[i]], getSumWindow(), getTelType( getTeltoAna()[i] ), i ); + } + else + { + copyIPRTelAveraged( getSumWindow(), getTeltoAna()[i], i ); + } + + cout << "*******************************" << endl; + cout << "* Starting IPR in Time Slices *" << endl; + cout << "*******************************" << endl; + if( true ) + { + cout << "MK ts: " << getCalData()->getPedsTS_vector( false ).size() << endl; + for( unsigned int ts = 0 ; ts < getCalData()->getPedsTS_vector( false ).size() ; ts++ ) + { + calculateIPRGraphsTimeSlices( fPedFileNameC[getTeltoAna()[i]], ts, getSumWindow(), getTelType( getTeltoAna()[i] ), i ); + } + } + } + return true; +} + +/* + * calculate IPR graphs and write them to disk + * + * (this is done per telescope type) + * + */ + +bool VIPRCalculator::copyIPRTelAveraged( unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel ) +{ + + TGraphErrors* i_IPRGraph = getIPRGraph( iSummationWindow, true ); + + setTelID( getTeltoAna()[0] ); + TGraphErrors* i_IPRGraph_Tel0 = getIPRGraph( iSummationWindow, true ); + + setTelID( getTeltoAna()[i_tel] ); + + if( !i_IPRGraph ) + { + cout << "VIPRCalculator::copyIPRTelAveraged: no IPR graph found for telescope type " << iTelType << endl; + return false; + } + if( !i_IPRGraph_Tel0 ) + { + cout << "VIPRCalculator::copyIPRTelAveraged: no IPR graph found for telescope type " << getTeltoAna()[0] << endl; + return false; + } + + for( Int_t i = 0; i < i_IPRGraph_Tel0->GetN() ; i++ ) + { + i_IPRGraph->SetPoint( i, i_IPRGraph_Tel0->GetPointX( i ), i_IPRGraph_Tel0->GetPointY( i ) ); + } + + i_IPRGraph->SetMinimum( 1 ); + i_IPRGraph->SetTitle( i_IPRGraph_Tel0->GetTitle() ); + i_IPRGraph->GetXaxis()->SetTitle( i_IPRGraph_Tel0->GetXaxis()->GetTitle() ); + i_IPRGraph->GetYaxis()->SetTitle( i_IPRGraph_Tel0->GetYaxis()->GetTitle() ); + i_IPRGraph->SetName( Form( "IPRcharge_TelType%d_SW%d", ( int )iTelType, iSummationWindow ) ); + + return true; + +} + +bool VIPRCalculator::calculateIPRGraphs( string iPedFileName, unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel ) +{ + cout << "MK sum window: " << iSummationWindow << endl; + TH1F* hIPR = ( TH1F* )initializeIPRHistogram( iSummationWindow, getTeltoAna()[i_tel] ); + TDirectory* iG_CurrentDirectory = gDirectory; + + // get an IPR graph + TGraphErrors* i_IPRGraph = getIPRGraph( iSummationWindow, true ); + if( !i_IPRGraph ) + { + cout << "VIPRCalculator::calculateIPRGraphs info: no IPR graph found for telescope type " << iTelType << endl; + return false; + } + // check suffix of ped file + if( iPedFileName.find( ".root" ) == string::npos ) + { + iPedFileName += ".root"; + } + // open pedestal files + TFile iPedFile( iPedFileName.c_str() ); + if( iPedFile.IsZombie() ) + { + cout << "VIPRCalculator::calculateIPRGraphs error reading IPR graphs from "; + cout << iPedFileName << endl; + return false; + } + // histograms with IPR distributions are either + // i) in a directory called distributions_TelType + // ii) in the current directory + cout << "Telescope type " << iTelType << ": "; + cout << "reading IPR histograms for summation window " << iSummationWindow; + cout << " from "; + cout << iPedFileName << endl; + + stringstream i_Directory( stringstream::in | stringstream::out ); + i_Directory << "distributions_" << iTelType; + if( iPedFile.Get( i_Directory.str().c_str() ) ) + { + iPedFile.cd( i_Directory.str().c_str() ); + } + + hIPR->Reset(); + + ////////////////////////////////////////////////// + // average over all channels in one telescope + float i_gainCorrect = 1.; + for( unsigned int i = 0; i < getNChannels(); i++ ) + { + if( getDetectorGeometry()->getAnaPixel()[i] > 0 + && i < getDead().size() && !getDead()[i] ) + { + stringstream i_Hname( stringstream::in | stringstream::out ); + i_Hname << "hpedPerTelescopeType_" << iTelType << "_" << iSummationWindow << "_" << i; + TH1F* h = ( TH1F* )gDirectory->Get( i_Hname.str().c_str() ); + + if( h ) + { + float ped = 0; + // default: pedestals are subtracted here + // (for combined channel analysis: charges are filled already pedestal subtracted) + // apply relative gain correction to integrated charges + if( getRunParameter()->fCombineChannelsForPedestalCalculation == 0 ) + { + ped = getPeds()[i]; + if(i<5){ cout << "MK ped: " << ped << endl; } + } + for( int j = 1; j <= h->GetNbinsX(); j++ ) + { + if( h->GetBinContent( j ) > 0. && getGains()[i] > 0 ) + { + i_gainCorrect = getGains()[i]; + if( getHIGHQE_gainfactor( i ) > 0. ) + { + i_gainCorrect *= getHIGHQE_gainfactor( i ); + } + if( i_gainCorrect > 0. ) + { + hIPR->Fill( ( h->GetBinCenter( j ) - ped * iSummationWindow ) / i_gainCorrect, h->GetBinContent( j ) ); + } + } + } + } + } + } + + int z = 0; + float norm = hIPR->Integral( 1, hIPR->GetNbinsX() ); + if( norm < fPedPerTelescopeTypeMinCnt ) //statistical limit for number of counts + { + fIPRAverageTel = true; + cout << "Telescope " << iTelType << ": "; + cout << "VIPRCalculator::calculateIPRGraphs WARNING: too few statistics to measure IPR curves "; + cout << "(total counts available: " << norm << ", "; + cout << "current limit " << fPedPerTelescopeTypeMinCnt << ")" << endl; + cout << "IPR graphs will be provided as sum of " << getTeltoAna().size() << " telescopes statistics." << endl; + cout << "VIPRCalculator::calculateIPRGraphs(): fIPRAverageTel = " << fIPRAverageTel << endl; + hIPR = calculateIPRGraphAveraged( iSummationWindow ); + cout << "VIPRCalculator::calculateIPRGraphs norm IPR combined: " << hIPR->Integral( 1, hIPR->GetNbinsX() ) << endl; + } + if( norm == 0 ) + { + cout << "VIPRCalculator::calculateIPRGraphs ERROR: no counts in IPR histogram !" << endl; + return false; + } + // convert to Rate + float nsToSec = 1E-9; + float Tsearch = getDetectorGeometry()->getLengthOfSampleTimeSlice( getTelID() ); + if( getSearchWindowLast() < getNSamples() ) + { + Tsearch *= ( getSearchWindowLast() - getSumFirst() ); // [ns] + } + else + { + Tsearch *= ( getNSamples() - getSumFirst() ); // [ns] + } + float convToHz = 1.; + if( nsToSec > 0. && Tsearch > 0. ) + { + convToHz /= ( nsToSec * Tsearch ); + } + else if( getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_ifExplicitSampleTimeSlice + && getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_sampleTimeSlice > 0 + && getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_nBinsADC > 0 ) + { + // simple peak sensing: sim_telarray uses the maximum bin only + // The values for sampleTimeSlice and nBinsADC are set in the cleaning parameter file + // For example, for the currect (Apr 17) ASTRI simulation, it is set in sim_telarray as + // fadc_mhz = 500 % MHz ==> sampleTimeSlice = 2 ns + // fadc_sum_bins = nBinsADC = 25 % Number of ADC time intervals actually summed up. + convToHz /= ( nsToSec + * getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_sampleTimeSlice + * getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_nBinsADC ); + } + + for( int i = 1; i <= hIPR->GetNbinsX(); i++ ) + { + if( hIPR->GetBinContent( i ) > 5 ) + { + double val = convToHz * hIPR->Integral( i, hIPR->GetNbinsX() ) / norm; + double valerr = convToHz * sqrt( hIPR->Integral( i, hIPR->GetNbinsX() ) ) / norm; + double charge_pe = hIPR->GetXaxis()->GetBinCenter( i ) * getTelescopeAverageFADCtoPhe(); + double charge_pe_bin_width = 0.5 * hIPR->GetXaxis()->GetBinWidth( i ) * getTelescopeAverageFADCtoPhe(); + + i_IPRGraph->SetPoint( z, charge_pe, val ); + i_IPRGraph->SetPointError( z, charge_pe_bin_width, valerr ); + z++; + } + } + + i_IPRGraph->SetMinimum( 1 ); + i_IPRGraph->SetTitle( Form( "Rate vs Threshold. W_{RO}=%2.1f ns, W_{int}=%2.1f ns", Tsearch, + getDetectorGeometry()->getLengthOfSampleTimeSlice( getTelID() ) * iSummationWindow ) ); + if( getRunParameter()->fIgnoreDSTGains ) + { + i_IPRGraph->GetXaxis()->SetTitle( "Threshold [FADC counts]" ); + } + else + { + i_IPRGraph->GetXaxis()->SetTitle( "Threshold [p.e.]" ); + } + i_IPRGraph->GetYaxis()->SetTitle( "Rate above Threshold [Hz]" ); + i_IPRGraph->SetName( Form( "IPRcharge_TelType%d_SW%d", ( int )iTelType, iSummationWindow ) ); + hIPR->Delete(); + + iPedFile.Close(); + + iG_CurrentDirectory->cd(); + + return true; +} + +TH1F* VIPRCalculator::initializeIPRHistogram( unsigned int iSummationWindow, unsigned int iTelType ) +{ + //get reference hpedPerTelescopeType for channel 1 + cout << "VIPRCalculator::initializeIPRHistogram initializing IPR graph for calculation of average of telescopes." << endl; + TH1F* hIPR = 0; + + + int TelType = getTelType( iTelType ); + string PedFile_name = fPedFileNameC[iTelType]; + + // check suffix of ped file + if( PedFile_name.find( ".root" ) == string::npos ) + { + PedFile_name += ".root"; + } + // open pedestal files + TFile PedFile( PedFile_name.c_str() ); + + if( PedFile.IsZombie() ) + { + cout << "VIPRCalculator::calculateIPRGraphsStatistics: error reading IPR graphs from "; + cout << PedFile_name << endl; + return hIPR; + } + + // histograms with IPR distributions are either + // i) in a directory called distributions_TelType + // ii) in the current directory + cout << "Telescope type " << TelType << ": "; + cout << "reading IPR histograms for summation window " << iSummationWindow; + cout << " from "; + cout << PedFile_name << endl; + + stringstream Directory( stringstream::in | stringstream::out ); + Directory << "distributions_" << TelType; + if( !PedFile.Get( Directory.str().c_str() ) ) + { + return hIPR; + } + + // get charge distribution for first channel as reference histogram + stringstream hIPRname( stringstream::in | stringstream::out ); + hIPRname << "distributions_" << TelType << "/hpedPerTelescopeType_" << TelType << "_" << iSummationWindow << "_" << 1; + + //? + TH1F* href = ( TH1F* )gDirectory->Get( hIPRname.str().c_str() ); + PedFile.Close(); + + if( !href ) + { + cout << " Error: cannot find IPR histogram " << hIPRname.str().c_str(); + cout << " in file:" << PedFile_name << " ... exiting " << endl; + return hIPR; + } + + + //summary histogram + if( getRunParameter()->fIgnoreDSTGains ) + { + // work in dc + hIPR = new TH1F( "", "", int( 1.5 * href->GetNbinsX() + 0.5 ), href->GetXaxis()->GetXmin(), href->GetXaxis()->GetXmax() ); + cout << "Error: " << href->GetNbinsX() << " " << href->GetXaxis()->GetXmin() << " " << href->GetXaxis()->GetXmax() << endl; + return hIPR; + } + else + { + // work in pe + hIPR = new TH1F( "", "", 1000, 0., 100. ); + return hIPR; + } +} + +TH1F* VIPRCalculator::calculateIPRGraphAveraged( unsigned int iSummationWindow ) +{ + + TH1F* hIPR = ( TH1F* )initializeIPRHistogram( iSummationWindow, getTeltoAna()[0] ); + + for( unsigned int teltype = 0; teltype < getTeltoAna().size(); teltype++ ) + { + setTelID( getTeltoAna()[teltype] ); + + // first find dead channels + // (ignore low gains here) + findDeadChans( false ); + // calculate IPR graphs + string PedFile_name = fPedFileNameC[getTeltoAna()[teltype]]; + int TelType = getTelType( getTeltoAna()[teltype] ); + + if( PedFile_name.find( ".root" ) == string::npos ) + { + PedFile_name += ".root"; + } + // open pedestal files + TFile PedFile( PedFile_name.c_str() ); + if( PedFile.IsZombie() ) + { + cout << "VIPRCalculator::calculateIPRGraphAveraged error reading IPR graphs from "; + cout << PedFile_name << endl; + TH1F* hNull = 0; + return hNull; + } + // histograms with IPR distributions are either + // i) in a directory called distributions_TelType + // ii) in the current directory + cout << "VIPRCalculator::calculateIPRGraphAveraged Telescope type " << TelType << ": "; + cout << "reading IPR histograms for summation window " << iSummationWindow; + cout << " from "; + cout << PedFile_name << endl; + stringstream Directory( stringstream::in | stringstream::out ); + Directory.str( std::string() ); + Directory << "distributions_" << TelType; + cout << Directory.str().c_str() << endl; + if( PedFile.Get( Directory.str().c_str() ) ) + { + PedFile.cd( Directory.str().c_str() ); + } + else + { + cout << "VIPRCalculator::calculateIPRGraphAveraged no directory: " << Directory.str().c_str() << endl; + } + + /////////////////////////// + // average over all channels in one telescope + float i_gainCorrect = 1.; + for( unsigned int i = 0; i < getNChannels(); i++ ) + { + if( getDetectorGeometry()->getAnaPixel()[i] > 0 + && i < getDead().size() && !getDead()[i] ) + { + stringstream HistName( stringstream::in | stringstream::out ); + HistName << "hpedPerTelescopeType_" << TelType << "_" << iSummationWindow << "_" << i; + TH1F* h = ( TH1F* )gDirectory->Get( HistName.str().c_str() ); + + + if( h ) + { + float ped = 0; + // default: pedestals are subtracted here + // (for combined channel analysis: charges are filled already pedestal subtracted) + // apply relative gain correction to integrated charges + if( getRunParameter()->fCombineChannelsForPedestalCalculation == 0 ) + { + ped = getPeds()[i]; + } + // special treatment for ASTRI telescopes + else if( getRunParameter()->fCombineChannelsForPedestalCalculation == 2 ) + { + stringstream Pname( stringstream::in | stringstream::out ); + Pname << "hped_" << TelType << "_" << iSummationWindow << "_" << i; + TH1F* hP = ( TH1F* )gDirectory->Get( Pname.str().c_str() ); + if( hP ) + { + ped = hP->GetMean(); + } + } + for( int j = 1; j <= h->GetNbinsX(); j++ ) + { + if( h->GetBinContent( j ) > 0. && getGains()[i] > 0 ) + { + i_gainCorrect = getGains()[i]; + if( getHIGHQE_gainfactor( i ) > 0. ) + { + i_gainCorrect *= getHIGHQE_gainfactor( i ); + } + if( i_gainCorrect > 0. ) + { + hIPR->Fill( ( h->GetBinCenter( j ) - ped * iSummationWindow ) / i_gainCorrect, h->GetBinContent( j ) ); + } + } + } + } + } + } + PedFile.Close(); + + } + hIPR->Scale( 1. / getTeltoAna().size() ); + float norm = getTeltoAna().size() * hIPR->Integral( 1, hIPR->GetNbinsX() ); + cout << "VIPRCalculator::calculateIPRGraphAveraged normalization of average IPR histogram " << norm; + cout << ". Returning IPR histogram." << endl; + if( norm < fPedPerTelescopeTypeMinCnt ) //statistical limit for number of counts + { + cout << "VIPRCalculator::calculateIPRGraphAveraged WARNING: there is NOT enough statistics "; + cout << " ( < " << fPedPerTelescopeTypeMinCnt << ") even when averaging over all telescopes." << endl; + return hIPR; + } + else + { + cout << "VIPRCalculator::calculateIPRGraphsStatistics: there is enough statistics to average over telescopes. " << endl; + return hIPR; + } +} + + +/* + * + * TIME SLICES CASE: + * + */ + +bool VIPRCalculator::clearHistos() +{ + fpedcal_histo_storage.clear(); + if( fpedcal_histo_storage.empty() ) + { + cout << "MK empty 1" << endl; + return false; + } + for( int telID = 0; telID < fpedcal_histo_storage.size(); telID++ ) + { + for( int i = 0; i < fpedcal_histo_storage[telID].size(); i++ ) + { + for( int j = 0; j < fpedcal_histo_storage[telID][i].size(); j++ ) + { + for( int k = 0; k < fpedcal_histo_storage[telID][i][j].size(); k++ ) + { + delete fpedcal_histo_storage[telID][i][j][k]; + } + fpedcal_histo_storage[telID][i][j].clear(); + } + fpedcal_histo_storage[telID][i].clear(); + } + fpedcal_histo_storage[telID].clear(); + } + return true; +} + +vector>>> VIPRCalculator::getStorageHist() +{ + return fpedcal_histo_storage; +} + +TH1F* VIPRCalculator::getIPRPedestalHisto( const int telID, const int ts, const int pixel, const int sw ) +{ + if( fpedcal_histo_storage.empty() ) + { + cout << "*** MK empty" << endl; + return nullptr; + } + if( fpedcal_histo_storage[telID][ts][pixel][sw] ) + { + return fpedcal_histo_storage[telID][ts][pixel][sw]; + } + else + { + return nullptr; + } +} + +void VIPRCalculator::fillIPRPedestalHisto( const int telID, const int NTimeSlices, const vector>>& fpedcal_histo ) +{ + //fpedcal_histo_storage[telID].push_back(fpedcal_histo[telID]); + if( true ) + { + vector< vector< TH1F* > > v1; + for( int p = 0 ; p < getNChannels(); p++ ) + { + vector< TH1F* > v2; + for( int sw = 0 ; sw < 3; sw++ ) + { + if( fpedcal_histo[telID][p][sw] ) + { + v2.push_back( ( TH1F* )fpedcal_histo[telID][p][sw]->Clone() ) ; + } + } + v1.push_back( v2 ); + v2.clear(); + } + fpedcal_histo_storage[telID].push_back( v1 ); + v1.clear(); + } +} + +bool VIPRCalculator::calculateIPRGraphsTimeSlices( string iPedFileName, int TS, unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel ) +{ + TH1F* hIPR = ( TH1F* )initializeIPRHistogram( iSummationWindow, getTeltoAna()[i_tel] ); + TDirectory* iG_CurrentDirectory = gDirectory; + + // get an IPR graph + TGraphErrors* i_IPRGraph = getIPRGraphTimeSlice( true, TS ); + if( !i_IPRGraph ) + { + cout << "VIPRCalculator::calculateIPRGraphsTimeSlices info: no IPR graph found for telescope type " << iTelType << endl; + return false; + } + // check suffix of ped file + if( iPedFileName.find( ".root" ) == string::npos ) + { + iPedFileName += ".root"; + } + // open pedestal files + TFile iPedFile( iPedFileName.c_str() ); + if( iPedFile.IsZombie() ) + { + cout << "VIPRCalculator::calculateIPRGraphsTimeSlices error reading IPR graphs from "; + cout << iPedFileName << endl; + return false; + } + // histograms with IPR distributions are either + // i) in a directory called distributions_TelType + // ii) in the current directory + cout << "Telescope type " << iTelType << ": "; + cout << "reading IPR histograms for summation window " << iSummationWindow; + cout << " from "; + cout << iPedFileName << endl; + + stringstream i_Directory( stringstream::in | stringstream::out ); + i_Directory << "distributions_" << iTelType; + if( iPedFile.Get( i_Directory.str().c_str() ) ) + { + iPedFile.cd( i_Directory.str().c_str() ); + } + + hIPR->Reset(); + + ////////////////////////////////////////////////// + // average over all channels in one telescope + float i_gainCorrect = 1.; + for( unsigned int i = 0; i < getNChannels(); i++ ) + { + if( getDetectorGeometry()->getAnaPixel()[i] > 0 + && i < getDead().size() && !getDead()[i] ) + { + stringstream i_Hname( stringstream::in | stringstream::out ); + i_Hname << "hpedTimeSlices_Tel" << iTelType << "_TS" << TS << "_Pix" << i << "_SW" << iSummationWindow; + TH1F* h = ( TH1F* )gDirectory->Get( i_Hname.str().c_str() ); + + if( h ) + { + float ped = 0; + // default: pedestals are subtracted here + // (for combined channel analysis: charges are filled already pedestal subtracted) + // apply relative gain correction to integrated charges + if( getRunParameter()->fCombineChannelsForPedestalCalculation == 0 ) + { + ped = getCalData()->getPedsTS_vector( false )[TS][i]; + if(i<5){ cout << "MK ped ts: " << ped << endl; } + } + for( int j = 1; j <= h->GetNbinsX(); j++ ) + { + if( h->GetBinContent( j ) > 0. && getGains()[i] > 0 ) + { + i_gainCorrect = getGains()[i]; + if( getHIGHQE_gainfactor( i ) > 0. ) + { + i_gainCorrect *= getHIGHQE_gainfactor( i ); + } + if( i_gainCorrect > 0. ) + { + hIPR->Fill( ( h->GetBinCenter( j ) - ped*iSummationWindow ) / i_gainCorrect, h->GetBinContent( j ) ); + } + } + } + } + } + } + + int z = 0; + float norm = hIPR->Integral( 1, hIPR->GetNbinsX() ); + cout << "MK norm of IPR: " << norm << endl; + if( norm < fPedPerTelescopeTypeMinCnt ) //statistical limit for number of counts + { + fIPRAverageTel = true; + cout << "Telescope " << iTelType << ": "; + cout << "VIPRCalculator::calculateIPRGraphsTimeSlices WARNING: too few statistics to measure IPR curves "; + cout << "(total counts available: " << norm << ", "; + cout << "current limit " << fPedPerTelescopeTypeMinCnt << ")" << endl; + cout << "IPR graphs will be provided as sum of " << getTeltoAna().size() << " telescopes statistics." << endl; + cout << "VIPRCalculator::calculateIPRGraphsTimeSlices(): fIPRAverageTel = " << fIPRAverageTel << endl; + //hIPR = calculateIPRGraphAveraged( iSummationWindow ); + cout << "VIPRCalculator::calculateIPRGraphsTimeSlices norm IPR combined: " << hIPR->Integral( 1, hIPR->GetNbinsX() ) << endl; + } + if( norm == 0 ) + { + cout << "VIPRCalculator::calculateIPRGraphsTimeSlices ERROR: no counts in IPR histogram !" << endl; + return false; + } + // convert to Rate + float nsToSec = 1E-9; + float Tsearch = getDetectorGeometry()->getLengthOfSampleTimeSlice( getTelID() ); + if( getSearchWindowLast() < getNSamples() ) + { + Tsearch *= ( getSearchWindowLast() - getSumFirst() ); // [ns] + } + else + { + Tsearch *= ( getNSamples() - getSumFirst() ); // [ns] + } + float convToHz = 1.; + if( nsToSec > 0. && Tsearch > 0. ) + { + convToHz /= ( nsToSec * Tsearch ); + } + else if( getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_ifExplicitSampleTimeSlice + && getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_sampleTimeSlice > 0 + && getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_nBinsADC > 0 ) + { + // simple peak sensing: sim_telarray uses the maximum bin only + // The values for sampleTimeSlice and nBinsADC are set in the cleaning parameter file + // For example, for the currect (Apr 17) ASTRI simulation, it is set in sim_telarray as + // fadc_mhz = 500 % MHz ==> sampleTimeSlice = 2 ns + // fadc_sum_bins = nBinsADC = 25 % Number of ADC time intervals actually summed up. + convToHz /= ( nsToSec + * getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_sampleTimeSlice + * getRunParameter()->fImageCleaningParameters[i_tel]->fNNOpt_nBinsADC ); + } + + for( int i = 1; i <= hIPR->GetNbinsX(); i++ ) + { + if( hIPR->GetBinContent( i ) > 5 ) + { + double val = convToHz * hIPR->Integral( i, hIPR->GetNbinsX() ) / norm; + double valerr = convToHz * sqrt( hIPR->Integral( i, hIPR->GetNbinsX() ) ) / norm; + double charge_pe = hIPR->GetXaxis()->GetBinCenter( i ) * getTelescopeAverageFADCtoPhe(); + double charge_pe_bin_width = 0.5 * hIPR->GetXaxis()->GetBinWidth( i ) * getTelescopeAverageFADCtoPhe(); + + i_IPRGraph->SetPoint( z, charge_pe, val ); + i_IPRGraph->SetPointError( z, charge_pe_bin_width, valerr ); + z++; + } + } + + i_IPRGraph->SetMinimum( 1 ); + i_IPRGraph->SetTitle( Form( "Rate vs Threshold. W_{RO}=%2.1f ns, W_{int}=%2.1f ns", Tsearch, + getDetectorGeometry()->getLengthOfSampleTimeSlice( getTelID() ) * iSummationWindow ) ); + if( getRunParameter()->fIgnoreDSTGains ) + { + i_IPRGraph->GetXaxis()->SetTitle( "Threshold [FADC counts]" ); + } + else + { + i_IPRGraph->GetXaxis()->SetTitle( "Threshold [p.e.]" ); + } + i_IPRGraph->GetYaxis()->SetTitle( "Rate above Threshold [Hz]" ); + i_IPRGraph->SetName( Form( "IPRcharge_TelType%d_TS%d_SW%d", ( int )iTelType, TS, iSummationWindow ) ); + hIPR->Delete(); + + iPedFile.Close(); + + iG_CurrentDirectory->cd(); + + return true; +} + +/* + * + * write IPR graphs to disk (one per telescope type) + * + */ +bool VIPRCalculator::writeIPRgraphs( map>> &hped_vec, string iFile ) +{ + TFile* fgraphs = 0; + if( iFile.size() == 0 ) + { + fgraphs = new TFile( getRunParameter()->fIPRdatabaseFile, "RECREATE" ); + } + else + { + fgraphs = new TFile( iFile.c_str(), "UPDATE" ); + } + if( fgraphs->IsZombie() ) + { + cout << "VCalibrator::writeIPRgraphs error opening IPR output file: " << endl; + cout << "\t" << fgraphs->GetName() << endl; + return false; + } + + // tree with conditions for these IRPs + TTree* header = new TTree( "IPRheader", "IPR parameters" ); + unsigned int sumwin = 0; // [FADC slices] + unsigned int Nsamples = 0; // [FADC slices] + float ROwin = 0.; // [ns] + float FADCslice = 0.; // [ns] + unsigned int SignalExtractor = 0; // [ Selected Extractor ] + header->Branch( "SignalExtractor", &SignalExtractor, "SignalExtractor/i" ); + header->Branch( "SummationWindow", &sumwin, "SummationWindow/i" ); + header->Branch( "Nsamples", &Nsamples, "Nsamples/i" ); + header->Branch( "FADCtimeSlice", &FADCslice, "FADCtimeSlice/F" ); + header->Branch( "ReadoutWindow", &ROwin, "ReadoutWindow/F" ); + + // one graph per telescope ID + map< ULong64_t, bool > iTelDone; + for( unsigned int i = 0; i < getDetectorGeometry()->getTelType_list().size(); i++ ) + { + iTelDone[getDetectorGeometry()->getTelType_list()[i]] = false; + } + + for( unsigned int i = 0; i < getNTel(); i++ ) + { + setTelID( i ); + if( iTelDone[getTelType( i )] ) + { + continue; + } + + if( hped_vec.find( getTelType( i ) ) == hped_vec.end() ) + { + continue; + } + + // loop over all summation windows + for( unsigned int j = 0; j < hped_vec[getTelType( i )].size(); j++ ) + { + // summation window + int i_sw = j + 1; + + TGraphErrors* g = getIPRGraph( i_sw, false ); + if( !g ) + { + continue; + } + SignalExtractor = getRunParameter()->fTraceIntegrationMethod.at( getTelID() ); + sumwin = i_sw; + Nsamples = getNSamples(); + FADCslice = getDetectorGeometry()->getLengthOfSampleTimeSlice( getTelID() ); //[ns] + ROwin = getDetectorGeometry()->getLengthOfSampleTimeSlice( getTelID() ) * ( float )getNSamples(); // [ns] + + header->Fill(); + g->Write(); + } + + //loop over all time slices + for (unsigned int ts = 0 ; ts < getCalData()->getPedsTS_vector( false ).size() ; ts++ ) + { + TGraphErrors* gTS = getIPRGraphTimeSlice(false, ts ); + if( !gTS ) + { + continue; + } + SignalExtractor = getRunParameter()->fTraceIntegrationMethod.at( getTelID() ); + sumwin = 1; + Nsamples = getNSamples(); + FADCslice = getDetectorGeometry()->getLengthOfSampleTimeSlice( getTelID() ); //[ns] + ROwin = getDetectorGeometry()->getLengthOfSampleTimeSlice( getTelID() ) * ( float )getNSamples(); // [ns] + + header->Fill(); + gTS->Write(); + + } + + iTelDone[getTelType( i )] = true; + } + header->Write(); + fgraphs->Close(); + fgraphs->Delete(); + + return true; +} diff --git a/src/VImageCleaning.cpp b/src/VImageCleaning.cpp index 27d52a18a..23bc72571 100644 --- a/src/VImageCleaning.cpp +++ b/src/VImageCleaning.cpp @@ -548,6 +548,19 @@ bool VImageCleaning::InitNNImgClnPerTelType( unsigned int teltype ) } // IPR graphs IPRgraph->Write(); + + for (unsigned int ts = 0 ; ts < 4 ; ts++ ) + { + TGraphErrors* gTS = fData->getIPRGraphTimeSlice(false, ts ); + if( !gTS ) + { + continue; + } + cout << "MK writing IPR" << endl; + gTS->Write(); + + } + if( fWriteGraphToFileRecreate ) { // probability curves (note that the x-axis is not charge!) diff --git a/src/VPedestalCalculator.cpp b/src/VPedestalCalculator.cpp index 1bafb52a2..4adc0411e 100644 --- a/src/VPedestalCalculator.cpp +++ b/src/VPedestalCalculator.cpp @@ -12,7 +12,7 @@ VPedestalCalculator::VPedestalCalculator() fDebug = getDebugFlag(); // default parameters (can be adjusted later in initialize()) - fLengthofTimeSlice = 180.; // in [s] + fLengthofTimeSlice = 2000.; // in [s] fSumWindow = 24; fNPixel = 500; fSumFirst = 0; @@ -184,9 +184,11 @@ bool VPedestalCalculator::initialize( bool ibCalibrationRun, unsigned int iNPixe fpedcal_mean.push_back( iped_cal2 ); fpedcal_mean2.push_back( iped_cal2 ); fpedcal_histo.push_back( iped_histo2 ); + fpedcal_histo_slidingw.push_back( iped_histo2 ); // define the time vector fTimeVec.push_back( 0 ); + NTimeSlices.push_back( 0 ); } iDir->cd(); @@ -252,6 +254,7 @@ void VPedestalCalculator::fillTimeSlice( unsigned int telID ) fpedcal_mean[telID][p][w] = 0.; fpedcal_mean2[telID][p][w] = 0.; fpedcal_histo[telID][p][w]->Reset(); + fpedcal_histo_slidingw[telID][p][w]->Reset(); } // deroate the pixel coordinates if( getTelID() < getPointing().size() && getPointing()[getTelID()] ) @@ -277,7 +280,7 @@ void VPedestalCalculator::fillTimeSlice( unsigned int telID ) } -void VPedestalCalculator::doAnalysis( bool iLowGain ) +void VPedestalCalculator::doAnalysis( bool iLowGain, VIPRCalculator *fIPRCalculator ) { double t = getEventTime(); // get right index for tel id @@ -304,12 +307,18 @@ void VPedestalCalculator::doAnalysis( bool iLowGain ) else if( t - fTimeVec[telID] > fLengthofTimeSlice ) { time = t; + + NTimeSlices[telID]+=1; + cout << "MK filling pedestal for IPR" << endl; + fIPRCalculator->fillIPRPedestalHisto(telID, NTimeSlices[telID], fpedcal_histo_slidingw ); + fillTimeSlice( telID ); fTimeVec[telID] = t; } // if( t - fTimeVec[telID] > fLengthofTimeSlice ) /////////////////////////////////////////////////////// double i_tr_sum = 0.; + double i_tr_sum_slidingw = 0.; // calculate the sums (don't use calcsums because it overwrites getSums() ) // and fill the histograms for( unsigned int i = 0; i < getNChannels(); i++ ) @@ -345,6 +354,7 @@ void VPedestalCalculator::doAnalysis( bool iLowGain ) { // calculate trace sum i_tr_sum = fTraceHandler->getTraceSum( fSumFirst, fSumFirst + ( w + 1 ), true, 1 ); + i_tr_sum_slidingw = fTraceHandler->getTraceSum( fSumFirst, fSumFirst + ( w + 1 ), false, 2 ); if( i_tr_sum > 0. && i_tr_sum < 50.*( w + 1 ) ) { if( chanID < fpedcal_n[telID].size() && w < fpedcal_n[telID][chanID].size() ) @@ -353,6 +363,7 @@ void VPedestalCalculator::doAnalysis( bool iLowGain ) fpedcal_mean[telID][chanID][w] += i_tr_sum; fpedcal_mean2[telID][chanID][w] += i_tr_sum * i_tr_sum; fpedcal_histo[telID][chanID][w]->Fill( i_tr_sum ); + fpedcal_histo_slidingw[telID][chanID][w]->Fill( i_tr_sum_slidingw ); } else { @@ -381,7 +392,7 @@ void VPedestalCalculator::doAnalysis( bool iLowGain ) } -void VPedestalCalculator::terminate( bool iWrite, bool iDebug_IO ) +void VPedestalCalculator::terminate( bool iWrite, bool iDebug_IO, VIPRCalculator* fIPRCalculator) { if( iWrite ) {