From 2961f0b048ec1f4043807aa23a5a2f50d2d4b2de Mon Sep 17 00:00:00 2001
From: Samir Amar-Youcef <samir@jspc54.x-matter.uni-frankfurt.de>
Date: Thu, 14 Nov 2013 13:32:30 +0100
Subject: [PATCH] Added: dynamic noise calculation in MAPS

---
 newCOMBI/MAPS.C | 232 +++++++++++++++++++++++++++++++++++++++++++++++-
 newCOMBI/MAPS.h |  13 ++-
 newCOMBI/run.C  |   7 +-
 3 files changed, 246 insertions(+), 6 deletions(-)

diff --git a/newCOMBI/MAPS.C b/newCOMBI/MAPS.C
index c95280c..bc53b67 100755
--- a/newCOMBI/MAPS.C
+++ b/newCOMBI/MAPS.C
@@ -87,10 +87,12 @@ MAPS::~MAPS(void) {
 		 
 		fHitTree->Write("",TObject::kOverwrite);
 		fNoiseTree->Write("",TObject::kOverwrite);
+		fDynNoiseTree->Write("",TObject::kOverwrite);
 		
 		hint1->Write("",TObject::kOverwrite);
 		hint2->Write("",TObject::kOverwrite);
 		
+		
 		fOutputFile->Save();
 		fOutputFile->Close();
 		
@@ -122,6 +124,8 @@ MAPS::~MAPS(void) {
 	delete[] fCdsmatrix;
 	delete[] fNoise;
 	delete[] fPedestals;
+	delete[] fDynFrameArr;
+	delete[] fDynCounter;
 	cout<<"================================================================="<<endl;
 	
 };
@@ -145,6 +149,7 @@ void MAPS::initMapsRun( TString InDir, TString OutDir, Int_t RunNumber, Int_t Ro
 	fOk			= false;
 	fNoiseOk	= false;
 	fNoiseExt	= false;
+	fNoiseDyn	= false;
 //-----------------------------------------------
 	cout<<"================================================================="<<endl;
 	cout<<"-----------------------"<<endl;
@@ -248,6 +253,7 @@ bool MAPS::checkConf( UInt_t &FileEvNb, UInt_t &TotEvNb, Int_t &PixelData ) {
 		fCdsmatrix	= new Int_t		[fPixelsData];
 		fNoise		= new Float_t	[fPixelsData];
 		fPedestals	= new Float_t	[fPixelsData];
+		fHittedPixel= new Int_t		[fPixelsData];
 //-----------------------------------------------
 		return true;
 	}
@@ -392,6 +398,7 @@ void MAPS::initSave() {
 		for(int i=0;i<25;i++)	{ fHitTree->Branch( Form("p%i",i+1)	, &fFrameInfo.p[i][0]	, Form("p%i [hits]/F",i+1)	, 32000); }
 // Noise TTree
 		fNoiseTree		= new TTree("noise", "noise");
+		fDynNoiseTree	= new TTree("dynnoise", "dynnoise");
 //-----------------------------------------------
 		hint1		= new TH2F("Hitmulit", "Hit multiplicity", 			fColumns, 0, fColumns, fRows, 0, fRows);
 		hint2		= new TH2F("Pixmulit", "Pix multiplicity", 			fColumns, 0, fColumns, fRows, 0, fRows);
@@ -554,7 +561,6 @@ bool MAPS::getNoise(Int_t Start, Int_t Frames) {
 			
 			for(Int_t j=0; j<Frames;j++)
 			{
-				if(ARR[j*fPixels+i]>5*PEDESTAL) continue; //throw all frames with a signal >5pedestal away, because of postulated signal charge
 				NOISE+=TMath::Power(ARR[j*fPixels+i]-PEDESTAL,2);
 			}
 
@@ -581,6 +587,193 @@ bool MAPS::getNoise(Int_t Start, Int_t Frames) {
 
 //####################################################################
 
+bool MAPS::getDynNoise(Int_t Frames) {
+//Dynamic Noise + Pedestals calculation
+	
+	if( fNoiseDyn )
+	{
+		delete[] fDynFrameArr;
+		delete[] fDynCounter;
+		
+		fDynFrameArr	= new Float_t[Frames*fPixels];
+		fDynCounter		= new Int_t[fPixels];
+	}
+	else
+	{
+		fDynNoiseTree->Branch("frame"		, &fFrameNumber		, "frame/i"		, 32000);
+		fDynNoiseTree->Branch("noise"		, &fNoiseMean		, "noise/F"		, 32000);
+		fDynNoiseTree->Branch("pedestal"	, &fPedestalsMean	, "pedestal/F"	, 32000);
+	}
+	
+	if(fOk)
+	{
+//Check amount of Frames
+		if(Frames<50)
+		{
+			cout<<"Too few Frames! Change Frames from "<<Frames<<" to ";
+			Frames = 50;
+			cout<<Frames<<" !"<<endl;
+		}
+// Warning: if too many pixels the dynamic array overflow the ram
+		Int_t End = Frames*fPixels*4;
+		cout<<"Run getDynNoise() initiated ..."<<endl;
+		cout<<"---"<<endl;
+		cout<<"'fDynFrameArr' allocates ";
+		if		(End/(TMath::Power(2,40)) >= 1)	{ printf("%6.2f TB  ", (Float_t)(End/(TMath::Power(2,40))) ); }
+		else if	(End/(TMath::Power(2,30)) >= 1)	{ printf("%6.2f GB  ", (Float_t)(End/(TMath::Power(2,30))) ); }
+		else if	(End/(TMath::Power(2,20)) >= 1)	{ printf("%6.2f MB  ", (Float_t)(End/(TMath::Power(2,20))) ); }
+		else if	(End/(TMath::Power(2,10)) >= 1)	{ printf("%6.2f kB  ", (Float_t)(End/(TMath::Power(2,10))) ); }
+		else									{ printf("%6.2f B   ", (Float_t)(End)); }
+		cout<<"for dynamic noise calculation! "<<endl;
+		
+		fDynFrameArr	= new Float_t[Frames*fPixels];
+		fDynCounter		= new Int_t[fPixels];
+		
+		for(Int_t i=0; i<fPixels;i++)	{ fDynCounter[i]=0; }
+		
+		bool STOP	= false;
+		fNoiseDyn	= false;
+			
+//Check availabilty of Frames
+		if( Frames			> fEventsSum )	{ Frames = fEventsSum;			printf("Changed 'Number Frames' to: %u\n", Frames ); }
+// Get fPixels and allocate array
+		Float_t	PEDESTAL;
+		Float_t	NOISE;
+// 		Int_t	PIXEL=0;
+		
+		if(!fNoiseOk) {getNoise(0,Frames);}
+		
+		{
+			int i=0;
+			
+			while(!STOP)
+			{
+				getFrame(i);
+				hitana();
+				
+				STOP = true;
+				
+				for(Int_t j=0; j<fPixels;j++)
+				{
+					if( fHittedPixel[j]==0)
+					{
+						if( fDynCounter[j]<Frames )
+						{
+							fDynFrameArr[Frames*j+fDynCounter[j]] = fCdsmatrix[j];//-fPedestals[j];
+							fDynCounter[j]++;
+							
+							STOP = false;
+						}
+					}
+				}
+				i++;
+			}
+		}
+		
+		for(Int_t i=0; i<fPixels;i++)
+		{
+			PEDESTAL	= 0;
+			NOISE		= 0;
+			
+			for(Int_t j=0; j<Frames;j++)
+			{
+				PEDESTAL+=fDynFrameArr[Frames*i+j];
+			}
+
+			PEDESTAL = PEDESTAL/Frames;
+			
+			for(Int_t j=0; j<Frames;j++)
+			{
+				NOISE+=TMath::Power(fDynFrameArr[Frames*i+j]-PEDESTAL,2);
+			}
+
+			NOISE = TMath::Sqrt(NOISE/(Frames));
+			
+			fNoise[i]		= NOISE;
+			fPedestals[i]	= PEDESTAL;
+// 			PIXEL			= i;
+			
+			fDynCounter[i] = 0;
+		}
+		
+		cout<<"-----------------------"<<endl;
+		
+		fNoiseOk	= true;
+		fNoiseDyn	= Frames;
+		
+		if( fSave  )
+		{
+			fFrameNumber	= 0;
+			fNoiseMean		= 0;
+			fPedestalsMean	= 0;
+			
+			fDynNoiseTree->Fill();
+			
+						fFrameNumber	= 0;
+			fNoiseMean		= TMath::Mean((const int)fPixels, fNoise);
+			fPedestalsMean	= TMath::Mean((const int)fPixels, fPedestals);
+			
+			fDynNoiseTree->Fill();
+		}
+		return true;
+	}
+	else
+	{
+			return false;
+	}
+};
+
+//####################################################################
+
+bool MAPS::regetDynNoise() {
+	
+	if( fNoiseDyn )
+	{
+		Float_t	PEDESTAL;
+		Float_t	NOISE;
+		Int_t	Frames = fNoiseDyn;
+		
+		for(Int_t i=0; i<fPixels;i++)
+		{
+			PEDESTAL	= 0;
+			NOISE		= 0;
+			
+			for(Int_t j=0; j<Frames;j++)
+			{
+				PEDESTAL+=fDynFrameArr[Frames*i+j];
+			}
+
+			PEDESTAL = PEDESTAL/Frames;
+			
+			for(Int_t j=0; j<Frames;j++)
+			{
+				NOISE+=TMath::Power(fDynFrameArr[Frames*i+j]-PEDESTAL,2);
+			}
+
+			NOISE = TMath::Sqrt(NOISE/(Frames));
+			
+			fNoise[i]		= NOISE;
+			fPedestals[i]	= PEDESTAL;
+		}
+		
+		if(fSave)
+		{	
+			fNoiseMean		= TMath::Mean((const int)fPixels, fNoise);
+			fPedestalsMean	= TMath::Mean((const int)fPixels, fPedestals);
+			
+			fDynNoiseTree->Fill();
+		}
+		
+		return true;
+	}
+	else
+	{
+		return false;
+	}
+}
+	
+//####################################################################
+
 bool MAPS::loadNoise( TString InDir, Int_t RunNumber, Int_t Matrix) {
 	
 	TString FILENAME = InDir+Form("%i/%i_%i.root",RunNumber,RunNumber,Matrix);
@@ -669,6 +862,8 @@ void MAPS::hitana() {
 	
 	for(Int_t i=0; i<fPixels; i++)
 	{
+		fHittedPixel[i] = 0;
+		
 		if( (float)(1.*fCdsmatrix[i]-fPedestals[i]) > (5.*fNoise[i]) )
 		{
 			HITNR++;
@@ -742,6 +937,26 @@ void MAPS::hitana() {
 			fHits++;
 			
 			if(fSave) { hint1->Fill( Hitlist[hit]%fColumns, (int)(Hitlist[hit]/fColumns) ); }
+			
+			for(Int_t row=0;row<5;row++)
+			{
+				for(Int_t column=0;column<5;column++)
+				{
+					if			( (row==0) && (A<2) )					{ }
+					else if 	( (row==1) && (A<1) )					{ }
+					else if 	( (row==3) && (A>= (fRows-1)) )			{ }
+					else if 	( (row==4) && (A>= (fRows-2)) )			{ }
+					else
+					{
+						if		( (column==0) && (B<2) )				{ }
+						else if	( (column==1) && (B==0) )				{ }
+						else if	( (column==3) && (B==fColumns-1))		{ }
+						else if	( (column==4) && (B>fColumns-3))		{ }
+						else											{ fHittedPixel[Hitlist[hit]+(row-2)*fColumns+(column-2)] = 1; }
+					}
+				}
+			}
+			fHittedPixel[Hitlist[hit]] = 2;
 		}
 	}
 //End: loop evaluate true seeds:
@@ -760,6 +975,21 @@ void MAPS::hitana() {
 		}
 	}
 	
+	if( fNoiseDyn )
+	{
+		for(Int_t j=0; j<fPixels;j++)
+		{
+			if( fHittedPixel[j]==0)
+			{
+				
+				fDynFrameArr[fNoiseDyn*j+fDynCounter[j]] = fCdsmatrix[j];
+				fDynCounter[j]++;
+				
+				if( fDynCounter[j] == fNoiseDyn ) { fDynCounter[j] = 0; }
+			}
+		}
+	}
+	
 	delete[] Hitlist;
 }}
 
diff --git a/newCOMBI/MAPS.h b/newCOMBI/MAPS.h
index c18f71c..8ddfd75 100644
--- a/newCOMBI/MAPS.h
+++ b/newCOMBI/MAPS.h
@@ -22,15 +22,21 @@ private:
 	bool		fFrameOk;
 	bool		fNoiseOk;
 	bool		fNoiseExt;
+	int			fNoiseDyn;
 	Float_t		HIT_tmp[29];
-		
+	
+	Float_t		fNoiseMean;
+	Float_t		fPedestalsMean;
 	Int_t*		fEvents;
 	Int_t*		fF0matrix;
 	Int_t*		fF1matrix;
 	Int_t*		fCdsmatrix;
 	Float_t*	fNoise;
 	Float_t*	fPedestals;
-		Int_t fHits;
+	Int_t fHits;
+	Float_t*	fDynFrameArr;
+	Int_t*		fDynCounter;
+	Int_t*		fHittedPixel;
 	
 	frameInfo fFrameInfo;// = new frameInfo;
 	
@@ -38,6 +44,7 @@ private:
 	TFile*		fOutputFile;
 	TTree*	fHitTree;
 	TTree*	fNoiseTree;
+	TTree*	fDynNoiseTree;
 	
 	TString fSystem;
 	
@@ -67,6 +74,8 @@ public:
 	
 	bool getFrame		(Int_t FrameNumber);
 	bool getNoise		(Int_t Start, Int_t Frames);
+	bool getDynNoise	(Int_t Frames);
+	bool regetDynNoise	();
 	bool loadNoise		(TString InDir, Int_t RunNumber, Int_t Matrix=0);
 	bool setNoise		(Float_t Noise);
 	bool setPedestals	(Float_t Pedestals);
diff --git a/newCOMBI/run.C b/newCOMBI/run.C
index 9703c87..1061874 100644
--- a/newCOMBI/run.C
+++ b/newCOMBI/run.C
@@ -27,7 +27,7 @@ void run()
 	Int_t	matrix		= 0;
 	bool	save		= 1;
 	
-	bool	loadNoise	= false;
+	bool	loadNoise		= false;
 	Int_t	runnumberNoise	= 29020;
 	
 //-------------------------------------------------------	
@@ -53,7 +53,8 @@ for(Int_t runLauf=0;runLauf<1;runLauf++)
 		}
 		else
 		{
-			a->getNoise		(0,1000);
+// 			a->getNoise		(0,1000);
+			a->getDynNoise	(1000);				// Dynamic Noise Calculation
 		}
 		
 		int start	= 0;
@@ -65,6 +66,7 @@ for(Int_t runLauf=0;runLauf<1;runLauf++)
 			a->getFrame(i);
 			a->filterCommonMode();
 			a->hitana();
+			if(i%100==0) { a->regetDynNoise(); }	// Dynamic Noise Calculation
 			
 			progress = (Int_t)(((i-start)*100)/(entries-1)*10);
 			if(progress!=progress_tmp)	{ print_progress( (((i-start)*100.)/(entries-1)) ); progress_tmp=progress;}
@@ -77,6 +79,5 @@ for(Int_t runLauf=0;runLauf<1;runLauf++)
 		
 		delete a;
 	//-------------------------------------------------------
-		cout<<entries<<endl;
 	}
 }
-- 
2.43.0