Time dependent IPR graphs for NN cleaning

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 \

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,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

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/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

inc/VIPRCalculator.h

@@ -0,0 +1,57 @@
+#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;
+		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:
+		vector<vector<vector<vector<TH1F*>>>> fpedcal_histo_storage;
+		bool calculateIPRGraphs( std::vector<std::string>  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<ULong64_t, vector<vector<TH1F*>>> &hped_vec, string iFile = "" );
+		void fillIPRPedestalHisto(const int telID, const int NTimeSlices,const vector<vector<vector<TH1F*>>>& fpedcal_histo );
+		void fillIPRPedestalHisto();
+		TH1F* getIPRPedestalHisto(const int telID, const int ts, const int pixel, const int sw);
+		bool clearHistos();
+		vector<vector<vector<vector<TH1F*>>>> getStorageHist();
+
+		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;
+		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,14 +79,14 @@ class VPedestalCalculator : public VImageBaseAnalyzer
 		VPedestalCalculator();
 		~VPedestalCalculator() {}
 
-		void doAnalysis( bool iLowGain = false );
+		void doAnalysis( bool iLowGain = false, VIPRCalculator *fIPRCalculator = 0);
 		vector< TTree* > getPedestalTree()
 		{
 			return fTree;
 		}
 		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

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,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++ )