Time dependent IPR graphs for NN cleaning

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 );

inc/VCalibrator.h

@@ -5,6 +5,11 @@
 
 #include "VImageBaseAnalyzer.h"
 #include "VPedestalCalculator.h"
+#include "VIPRCalculator.h"
+#include "VDB_CalibrationInfo.h"
+#include <VSQLTextFileReader.h>
+#include "VIPRCalculator.h"
+#include "VPedestalCalculator.h"
 #include "VDB_CalibrationInfo.h"
 #include "VSQLTextFileReader.h"
 
@@ -109,6 +114,7 @@ class VCalibrator : public VImageBaseAnalyzer
 		bool copyIPRTelAveraged( unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel );
 		TH1F* initializeIPRAveraged( unsigned int iSummationWindow, unsigned int iTelType );
 		TH1F* calculateIPRGraphAveraged( unsigned int iSummationWindow );
+		bool calculateIPRGraphs( string iPedFileName, unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel );
 		bool readLowGainMultiplier( );
 		bool readPeds( string iFile, bool, unsigned int );
 		bool readPeds_from_grisufile( bool, unsigned int );
@@ -122,6 +128,7 @@ class VCalibrator : public VImageBaseAnalyzer
 
 		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 +142,10 @@ class VCalibrator : public VImageBaseAnalyzer
 		void calculatePedestals( bool iLowGain = false );
 		void calculateGainsAndTOffsets( bool iLowGain = false );
 		unsigned int getNumberOfEventsUsedInCalibration( int iTelID, int iType );
-		void initialize();
+		void initialize(VIPRCalculator *fIPRCalculator );
 		void terminate( VPedestalCalculator* );
+		//void calculateIPRGraphs(VIPRCalculator *fIPRCalculator = 0, VPedestalCalculator *iP = 0);
+		void initialize( VIPRCalculator* fIPRCalculator );
+		void terminate( VPedestalCalculator* iP, VIPRCalculator* fIPRCalculator );
 };
 #endif

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

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 )

inc/VIPRCalculator.h

@@ -0,0 +1,73 @@
+#ifndef VIPRCALCULATOR_H
+#define VIPRCALCULATOR_H
+
+#include <VImageBaseAnalyzer.h>
+#include <VDB_CalibrationInfo.h>
+#include <VSQLTextFileReader.h>
+
+#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 <fstream>
+#include <map>
+#include <sstream>
+#include <string>
+#include <vector>
+
+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;
+<<<<<<< HEAD
+		TH1F* FillIPRHistogram( unsigned int iSummationWindow,  unsigned int i_tel);
+		void definePedestalFile( std::vector<std::string> 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:
+		bool calculateIPRGraphs( std::vector<std::string>  fPedFileNameCalibrator );
+		bool calculateIPRGraphs( string iPedFileName, unsigned int iSummationWindow, ULong64_t iTelType, unsigned int i_tel );
+		bool writeIPRgraphs( map<ULong64_t, vector<vector<TH1F*>>> &hped_vec, string iFile = "" );
+
+=======
+		TH1F* FillIPRHistogram( unsigned int iSummationWindow,  unsigned int i_tel );
+		TGraphErrors* updateIPRGraph( TH1F* hIPR, unsigned int i_tel, int iSummationWindow );
+		float convertRate( unsigned int i_tel );
+		float getTsearch();
+		void definePedestalFile( std::vector<std::string> fPedFileNameC );
+		TFile* initializePedestalFile( int i_tel );
+		TH1F* initializeIPRHistogram( unsigned int iSummationWindow, unsigned int i_tel );
+		bool calculateIPRGraphsTimeSlices( const  int TimeSlice, int iSummationWindow, const int itel );
+		bool calculateIPRGraphsNoTimeSlices( unsigned int iSummationWindow, unsigned int i_tel );
+		bool copyIPRInitialized( unsigned int iSummationWindow, unsigned int i_tel );
+		TH1F* calculateIPRGraphAveragedNoTimeSlices( unsigned int iSummationWindow, unsigned int i_tel );
+	public:
+
+		vector<vector<vector<vector<TH1F*>>>> fpedcal_histo_storage;
+		TH1F* getIPRPedestalHisto( const int telID, const int ts, const int pixel, const int sw );
+		vector<vector<vector<vector<TH1F*>>>> getStorageHist();
+		void fillIPRPedestalHisto();
+		void fillIPRPedestalHisto( const int telID, const vector<vector<vector<TH1F*>>>& fpedcal_histo );
+		bool calculateIPRGraphs( std::vector<std::string>  fPedFileNameC );
+		bool writeIPRgraphs( map<ULong64_t, vector<vector<TH1F*>>>& hped_vec, string iFile = "" );
+		void checkHistEmpty( const int telID, const int ts, const int pixel, const int sw );
+		VIPRCalculator();
+		~VIPRCalculator() {}
+
+		void initialize();
+};
+#endif

inc/VIPRCalculator.h.bsk

@@ -0,0 +1,62 @@
+#ifndef VIPRCALCULATOR_H
+#define VIPRCALCULATOR_H
+
+#include <VImageBaseAnalyzer.h>
+#include <VDB_CalibrationInfo.h>
+#include <VSQLTextFileReader.h>
+
+#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 <fstream>
+#include <map>
+#include <sstream>
+#include <string>
+#include <vector>
+
+using namespace std;
+
+class VIPRCalculator : public VImageBaseAnalyzer
+{
+
+	private:
+		VEvndispData*      fData;
+		vector< string > fPedFile;
+		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);
+		TGraphErrors* updateIPRGraph(TH1F *hIPR, unsigned int i_tel, int iSummationWindow);
+		float convertRate(unsigned int i_tel);
+		float getTsearch();
+		void definePedestalFile( std::vector<std::string> fPedFileNameC );
+		TFile* initializePedestalFile( int i_tel );
+		TH1F* initializeIPRHistogram( unsigned int iSummationWindow, unsigned int i_tel);
+		bool calculateIPRGraphsTimeSlices( const  int TimeSlice, int iSummationWindow, const int itel);
+		bool calculateIPRGraphsNoTimeSlices( unsigned int iSummationWindow, unsigned int i_tel );
+		bool copyIPRInitialized( unsigned int iSummationWindow, unsigned int i_tel );
+		TH1F* calculateIPRGraphAveragedNoTimeSlices( unsigned int iSummationWindow, unsigned int i_tel );
+	public:
+		void clearHistos();
+		vector<vector<vector<vector<TH1F*>>>> fpedcal_histo_storage;
+		TH1F* getIPRPedestalHisto(const int telID, const int ts, const int pixel, const int sw);
+		vector<vector<vector<vector<TH1F*>>>> getStorageHist();
+		void fillIPRPedestalHisto();
+		void fillIPRPedestalHisto(const int telID, const int NTimeSlices, const vector<vector<vector<TH1F*>>>& fpedcal_histo );
+		bool calculateIPRGraphs( std::vector<std::string>  fPedFileNameC);
+		bool writeIPRgraphs( map<ULong64_t, vector<vector<TH1F*>>> &hped_vec, string iFile = "" );
+		void checkHistEmpty(const int telID, const int ts, const int pixel, const int sw);
+
+		VIPRCalculator();
+                ~VIPRCalculator() {}
+
+		void initialize();
+};
+#endif 

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;
+		std::vector<std::vector<TH1F*> > copy_fpedcal_histo;
 
 		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]
@@ -78,6 +80,12 @@ class VPedestalCalculator : public VImageBaseAnalyzer
 		~VPedestalCalculator() {}
 
 		void doAnalysis( bool iLowGain = false );
+		vector< vector< vector< TH1F* > > > fpedcal_histo_sw;
+
+		VPedestalCalculator();
+		~VPedestalCalculator() {}
+
+		void doAnalysis( bool iLowGain = false , VIPRCalculator* fIPRCalculator = 0 );
 		vector< TTree* > getPedestalTree()
 		{
 			return fTree;

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;

src/VCalibrator.cpp

@@ -359,7 +359,9 @@ void VCalibrator::calculatePedestals( bool iLowGain )
    this might be an ascii file and/or a root file
 
 */
+//MK test
 void VCalibrator::writePeds( bool iLowGain, VPedestalCalculator* iPedestalCalculator, bool iWriteAsciiFile )
+//void VCalibrator::writePeds( bool iLowGain, VPedestalCalculator* iPedestalCalculator, bool iWriteAsciiFile, VIPRCalculator *fIPRCalculator )
 {
 	if( getDebugFlag() )
 	{
@@ -1588,12 +1590,14 @@ void VCalibrator::writeTOffsets( bool iLowGain )
 	}
 }
 
-
+//MK test
+//void VCalibrator::terminate( VPedestalCalculator* iP, VIPRCalculator *fIPRCalculator )
 void VCalibrator::terminate( VPedestalCalculator* iP )
 {
 	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 +3094,8 @@ void VCalibrator::readTOffsets( bool iLowGain )
 
 }
 
-
-void VCalibrator::initialize()
+//MK test
+void VCalibrator::initialize(VIPRCalculator *fIPRCalculator )
 {
 	if( fDebug )
 	{
@@ -3221,6 +3225,8 @@ void VCalibrator::initialize()
 	// if needed: write IPR graphs to disk
 	if( getRunParameter()->ifCreateIPRdatabase == true && getRunParameter()->ifReadIPRfromDatabase == false )
 	{
+		//MK test
+		//fIPRCalculator->writeIPRgraphs( hped_vec, "");
 		writeIPRgraphs();
 	}
 
@@ -3237,7 +3243,10 @@ void VCalibrator::initialize()
 				&& fRunPar->frunmode != 1
 				&& fRunPar->frunmode != 6 )
 		{
-			calculateIPRGraphs();
+			cout << "MK calculating IPR graphs"<< endl;
+			//calculateIPRGraphs();
+			fIPRCalculator->calculateIPRGraphs(fPedFileNameC);
+			//calculateIPRGraphs();
 		}
 	}
 }
@@ -4376,7 +4385,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++ )