Context Navigation

svv2mtx.cc@ 4059

Visit:

Last change on this file since 4059 was 3938, checked in by ansari, 15 years ago
modification provisoire pour lecture fichiers visibilites 2010, Reza 13/01/2011
File size: 13.6 KB

Rev	Line
[3646]	1	// Utilisation de SOPHYA pour faciliter les tests ...
	2	#include "sopnamsp.h"
	3	#include "machdefs.h"
	4
	5	/* ----------------------------------------------------------
[3701]	6	Projet BAORadio - (C) LAL/IRFU 2008-2010
	7
	8	Programme de lecture des fichiers vecteurs/matrices de
	9	visibilites produits par mcrd / vismfib
	10	R. Ansari, C. Magneville - LAL/Irfu
[3646]	11	V : Mai 2009
	12	---------------------------------------------------------- */
	13
	14	// include standard c/c++
	15	#include <math.h>
	16	#include <stdio.h>
	17	#include <stdlib.h>
	18	#include <string.h>
	19
	20	#include <iostream>
	21	#include <string>
	22
	23	#include "pexceptions.h"
	24	#include "tvector.h"
	25	#include "fioarr.h"
	26	// #include "tarrinit.h"
	27	#include "ntuple.h"
[3701]	28	#include "datatable.h"
[3646]	29	#include "histinit.h"
	30	#include "matharr.h"
	31	#include "timestamp.h"
[3701]	32	#include <utilarr.h>
[3646]	33
	34	// include sophya mesure ressource CPU/memoire ...
	35	#include "resusage.h"
	36	#include "ctimer.h"
	37	#include "timing.h"
	38
	39
[3701]	40	//--------------------------- Fonctions de ce fichier -------------------
[3646]	41	int Usage(bool fgshort=true);
	42
[3701]	43	int DecodeProc(int narg, char* arg[]);
	44	int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq,
	45	vector<sa_size_t> tfrlist);
[3704]	46	int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_4> chanum_,
[3701]	47	sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_) ;
	48
	49	int DecodeProcVJun09(int narg, char* arg[]);
	50	int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card);
	51	//------------------------------------------------------------------------------------------------------------
	52
	53	/* --Fonction-- */
	54	int Usage(bool fgshort)
	55	{
	56	cout << " --- svv2mtx.cc : Read PPF files produced by mcrd/visfmib to make time-frequency\n"
	57	<< " matrices and Visibilites=f(time) datatable " << endl;
	58	cout << " Usage: svv2mtx -mcrd InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq [card=1] \n"
	59	<< " OR \n"
	60	<< " svv2mtx InOutPath Imin,Imax,step NumFreq1,NumFreq2,NBinFreq VisMtxRowList" << endl;
	61	if (fgshort) {
	62	cout << " svv2mtx -h for detailed instructions" << endl;
	63	return 1;
	64	}
	65	cout << " -mcrd : Read mcrd output files -> ProcSVFilesVJun09() \n"
	66	<< " InOutPath : Input/Output directory name \n"
	67	<< " Imin,Imax,IStep: Input PPF files sequence number \n"
	68	<< " FileNames=InOutPath/Ch12_II.fits Imin<=II<=Imax II+=IStep \n"
	69	<< " NumFreq1,NumFreq2,NBinFreq: Freq Zone and number of frequency bins for ntuple\n"
	70	<< " VisMtxRowList : List of visibiliy matrix rows (0,2,...) -> time-freq map \n"
	71	<< " card=1 Ch12 , card=2 Ch34 (mcrd output files) " << endl;
	72	return 1;
	73	}
	74
[3646]	75	//----------------------------------------------------
	76	//----------------------------------------------------
	77	int main(int narg, char* arg[])
	78	{
	79	if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
	80	if (narg<4) return Usage(true);
	81	HiStatsInitiator _inia;
	82	// TArrayInitiator _inia;
	83
	84	int rc = 0;
	85	try {
[3701]	86	bool fgvjun09=false;
	87	if ((strcmp(arg[1],"-mcrd")==0)) fgvjun09=true;
[3646]	88	ResourceUsage resu;
[3701]	89	if (fgvjun09) DecodeProcVJun09(narg-2, arg+2);
	90	else DecodeProc(narg-1, arg+1);
[3646]	91	resu.Update();
	92	cout << resu;
	93	}
	94	catch (PException& exc) {
[3702]	95	cerr << " svv2mtx.cc catched PException " << exc.Msg() << endl;
[3646]	96	rc = 77;
	97	}
	98	catch (std::exception& sex) {
	99	cerr << "\n svv2mtx.cc std::exception :"
	100	<< (string)typeid(sex).name() << "\n msg= "
	101	<< sex.what() << endl;
	102	rc = 78;
	103	}
	104	catch (...) {
	105	cerr << " svv2mtx.cc catched unknown (...) exception " << endl;
	106	rc = 79;
	107	}
	108
	109	cout << ">>>> svv2mtx.cc ------- END ----------- RC=" << rc << endl;
	110	return rc;
	111
	112	}
	113
[3701]	114	//--------------------------------------------------------------------
	115	// Traitement fichiers produits par vismfib (V Nov09)
	116	/* --Fonction-- */
	117	int DecodeProc(int narg, char* arg[])
	118	{
	119	// Decodage des arguments et traitement
	120	string inoutpath = arg[0];
	121	int imin=0;
	122	int imax=0;
	123	int istep=1;
	124	sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
	125	int jf1=0;
	126	int jf2=0;
	127	int nfreq=0;
	128	sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
	129	int card=1;
	130	vector<sa_size_t> rowlist;
	131	if (narg>3) {
	132	EnumeratedSequence es;
	133	int nbad;
	134	es.Append(arg[3], nbad, ",");
	135	for(int j=0; j<es.Size(); j++) rowlist.push_back((int)es.Value(j));
	136	}
	137	if (rowlist.size()<1) rowlist.push_back(0);
[3646]	138
[3701]	139	cout << " ----- svv2mtx/DecodeProc - Start - InOutPath= " << inoutpath << " IMin,Max,Step="
	140	<< imin << "," << imax << "," << istep << " Card=" << card << endl;
	141	cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
	142	cout << " RowList: " ;
	143	for(int j=0; j<rowlist.size(); j++) cout << rowlist[j] << " , " ;
	144	cout << endl;
	145	int rc=ProcSVFiles(inoutpath, imin, imax, istep, jf1, jf2, nfreq, rowlist);
	146	return rc;
	147	}
[3646]	148
	149	// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
[3701]	150	/* --Fonction-- */
	151	int ProcSVFiles(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq,
	152	vector<sa_size_t> rowlist)
	153	{
	154	Timer tm("ProcSVFiles");
[3703]	155
[3701]	156	DataTable visdt;
	157	visdt.AddDoubleColumn("mfc");
	158	visdt.AddDoubleColumn("mtt");
	159	visdt.AddIntegerColumn("jfreq");
	160	visdt.AddIntegerColumn("numch");
	161	visdt.AddFloatColumn("vre");
	162	visdt.AddFloatColumn("vim");
	163
[3703]	164	vector< TMatrix< complex<r_4> > > vmtf;
	165
[3701]	166	if (jf1<1) jf1=1;
	167	if ((jf2<1)\|\|(jf2<jf1)) jf2=jf1;
	168	if (nfreq<1) nfreq=1;
	169	int djf=(jf2-jf1)/nfreq;
[3714]	170	if (djf<1) djf=1;
[3703]	171
[3701]	172	char fname[1024];
	173
[3703]	174	cout << " ---- ProcSVFiles()-Begin - Reading chanum vector" << endl;
[3938]	175	/*
	176	TVector< uint_4 > chanum(8*17);
	177	int ku=0;
	178	for(int iu=1; iu<=16; iu++)
	179	for(int ju=iu; ju<=16; ju++) {
	180	chanum(ku)=iu*1000+ju; ku++;
	181	}
	182	*/
	183	TVector< uint_4 > chanum;
[3703]	184	{
[3938]	185	sprintf(fname, "%s/chanum_0.ppf",inoutpath.c_str());
[3701]	186	PInPersist pic(fname);
[3938]	187	pic >> PPFNameTag("chanpairnum") >> chanum;
[3701]	188	}
[3938]	189
[3703]	190	cout << " ---- ProcSVFiles()-Read chanum done " << endl;
	191	sa_size_t nrows = (imax-imin+1)/istep;
[3701]	192	sa_size_t kr=0;
	193
[3938]	194	sa_size_t mtf_binfreq=25;
	195	sa_size_t mtf_bintime=5;
	196
	197	sa_size_t ncols=1;
	198
[3701]	199	for(int ifile=imin; ifile<=imax; ifile+=istep) {
[3938]	200	sprintf(fname, "%s/vismtx_0_%d.ppf",inoutpath.c_str(),ifile);
[3701]	201	cout << " ProcSVFiles[" << ifile << "] opening file " << fname << endl;
	202	PInPersist pin(fname);
	203	TMatrix< complex<r_4> > vismtx;
	204	pin >> vismtx;
	205
	206	if (ifile==imin) {
[3938]	207	ncols = vismtx.NCols();
[3701]	208	cout << " ProcSVFilesVJun09/Info: Output Time-Frequency matrices NRows=NFiles"
	209	<< nrows << " NCols=NFreq=" << ncols << endl;
[3703]	210	for(size_t j=0; j<rowlist.size(); j++)
[3938]	211	vmtf.push_back(TMatrix< complex<r_4> >(ncols/mtf_binfreq+1, nrows/mtf_bintime+1));
[3701]	212	}
[3938]	213	/*
[3701]	214	for(size_t j=0; j<rowlist.size(); j++)
	215	vmtf[j].Row(kr) = vismtx.Row(rowlist[j]);
[3938]	216	*/
	217	for(size_t j=0; j<rowlist.size(); j++)
	218	for(sa_size_t icf=0; icf<ncols; icf++) {
	219	vmtf[j](icf/mtf_binfreq,kr/mtf_bintime)+=vismtx(rowlist[j],icf);
	220	}
	221	// cout << " DBG* kr=" << kr << " kr/mtf_bintime=" << kr/mtf_bintime
	222	// << " ncols/2/mtf_binfreq=" << ncols/2/mtf_binfreq << endl;
[3701]	223	kr++;
[3938]	224
[3701]	225
	226	// Calcul moyenne dans des bandes en frequence
	227	FillVisDTable(visdt, vismtx, chanum, jf1, jf2, djf);
	228	}
	229
	230	sprintf(fname, "%s/vistfreqmtx.ppf",inoutpath.c_str());
	231	cout << "ProcSVFiles: Opening file " << fname << " for writing Visi(Freq,Time) matrices" << endl;
	232	POutPersist po(fname);
	233	char tagb[64];
	234	for(size_t j=0; j<rowlist.size(); j++) {
	235	sprintf(tagb,"visft%d", chanum(rowlist[j]));
	236	po << PPFNameTag(tagb) << vmtf[j];
	237	}
	238	cout << visdt;
	239	sprintf(fname, "%s/svvdt.ppf",inoutpath.c_str());
	240	cout << "ProcSVFile: writing visibility datatable to file " << fname << endl;
	241	POutPersist pod(fname);
	242	pod << visdt;
	243
	244	return 0;
	245	}
	246
	247	/* --Fonction-- */
	248	int FillVisDTable(BaseDataTable& visdt_, TMatrix< complex<r_4> > vismtx_, TVector< uint_4> chanum_,
	249	sa_size_t jf1_, sa_size_t jf2_, sa_size_t djf_)
	250	{
	251	double xnt_[20];
	252	xnt_[0]=(double)vismtx_.Info()["MeanFC"];
	253	xnt_[1]=(double)vismtx_.Info()["MeanTT"]/1.25e8;
	254
	255	uint_4 nmean_=vismtx_.Info()["NPAQSUM"];
	256	if (djf_<2) {
	257	for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
	258	for(sa_size_t jf=jf1_; jf<jf2_; jf++) {
	259	xnt_[2]=jf;
	260	xnt_[3]=chanum_(rv);
	261	xnt_[4]=vismtx_(rv,jf).real()/(r_4)(nmean_);
	262	xnt_[5]=vismtx_(rv,jf).imag()/(r_4)(nmean_);
	263	visdt_.AddRow(xnt_);
	264	}
	265	}
	266	}
	267	else {
	268	for(sa_size_t rv=0; rv<vismtx_.NRows(); rv++) {
	269	for(sa_size_t jf=jf1_; jf<jf2_; jf+=djf_) {
	270	r_4 moyreal=0.;
	271	r_4 moyimag=0.;
[3705]	272	sa_size_t jjfmx=jf+djf_;
	273	if (jjfmx > vismtx_.NCols()) jjfmx=vismtx_.NCols();
	274	for(sa_size_t jjf=jf; jjf<jjfmx; jjf++) {
[3701]	275	moyreal+=vismtx_(rv,jjf).real();
	276	moyimag+=vismtx_(rv,jjf).imag();
	277	}
	278	xnt_[2]=jf+djf_/2;
	279	xnt_[3]=chanum_(rv);
	280	xnt_[4]=moyreal/(r_4)(nmean_*djf_);
	281	xnt_[5]=moyimag/(r_4)(nmean_*djf_);
	282	visdt_.AddRow(xnt_);
	283	}
	284	}
	285	}
	286	return 0;
	287	}
	288
	289
	290	//--------------------------------------------------------------------
	291	// Traitement fichiers produits par mcrd (V Jun09)
	292	/* --Fonction-- */
	293	int DecodeProcVJun09(int narg, char* arg[])
	294	{
	295	// Decodage des arguments et traitement
	296	string inoutpath = arg[0];
	297	int imin=0;
	298	int imax=0;
	299	int istep=1;
	300	sscanf(arg[1],"%d,%d,%d",&imin,&imax,&istep);
	301	int jf1=0;
	302	int jf2=0;
	303	int nfreq=0;
	304	sscanf(arg[2],"%d,%d,%d",&jf1,&jf2,&nfreq);
	305	int card=1;
	306	if (narg>3) card=atoi(arg[3]);
	307	cout << " ----- svv2mtx/DecodeProcVJun09 - Start - InOutPath= " << inoutpath << " IMin,Max,Step="
	308	<< imin << "," << imax << "," << istep << " Card=" << card << endl;
	309	cout << "Frequency num range JF=" << jf1 << "," << jf2 << "," << nfreq << " ------- " << endl;
	310	int rc=ProcSVFilesVJun09(inoutpath, imin, imax, istep, jf1, jf2, nfreq, card);
	311	return rc;
	312	}
	313
	314
	315	// Pour traitement (calcul FFT et visibilites (ProcA) 1 fibre, 2 voies RAW)
	316	/* --Fonction-- */
[3699]	317	int ProcSVFilesVJun09(string& inoutpath, int imin, int imax, int istep, int jf1, int jf2, int nfreq, int card)
[3646]	318	{
	319	Timer tm("ProcSVFiles");
	320	char fname[512];
	321	// NTuple
[3647]	322	const char* nnames[10] = {"fcsm","ttsm","jfreq","s1","s2","s12","s12re","s12im","s12phi","s12mod"};
	323	NTuple nt(10, nnames);
	324	double xnt[15];
[3646]	325	uint_4 nmnt = 0;
[3647]	326	double ms1,ms2,ms12,ms12re,ms12im,ms12phi,ms12mod;
[3646]	327
	328	TMatrix<r_4> s1, s2;
[3647]	329	TMatrix<r_4> v12re, v12im, v12phi,v12mod;
[3646]	330	sa_size_t ncols = (imax-imin+1)/istep;
	331	sa_size_t nrows = 10;
	332	sa_size_t kc=0;
	333	for(int ifile=imin; ifile<=imax; ifile+=istep) {
	334	if (card==2)
	335	sprintf(fname, "%s/Ch34_%d.ppf",inoutpath.c_str(),ifile);
	336	else
	337	sprintf(fname, "%s/Ch12_%d.ppf",inoutpath.c_str(),ifile);
[3701]	338	cout << " ProcSVFilesVJun09[" << ifile << "] opening file " << fname << endl;
[3646]	339	PInPersist pin(fname);
	340	string tag1="specV1";
	341	string tag2="specV2";
	342	string tag12="visiV12";
	343	if (card==2) {
	344	tag1 = "specV3";
	345	tag2 = "specV4";
	346	tag12="visiV34";
	347	}
	348	TVector<r_4> sv1;
	349	TVector<r_4> sv2;
	350	TVector< complex<r_4> > vv12;
	351	pin >> PPFNameTag(tag1) >> sv1;
	352	pin >> PPFNameTag(tag2) >> sv2;
	353	pin >> PPFNameTag(tag12) >> vv12;
	354	if (ifile==imin) {
	355	nrows = sv1.Size();
[3699]	356	cout << " ProcSVFilesVJun09/Info: Output s1,s2 matrix size NRows=NFreq="
[3646]	357	<< nrows << " NCols=NFiles=" << ncols << endl;
	358	s1.SetSize(nrows, ncols);
	359	s2.SetSize(nrows, ncols);
	360	nrows = vv12.Size();
[3699]	361	cout << " ProcSVFilesVJun09/Info: Output v12 matrix size NRows=NFreq="
[3646]	362	<< nrows << " NCols=NFiles=" << ncols << endl;
	363	v12re.SetSize(nrows, ncols);
	364	v12im.SetSize(nrows, ncols);
	365	v12phi.SetSize(nrows, ncols);
[3647]	366	v12mod.SetSize(nrows, ncols);
[3646]	367	}
	368	s1.Column(kc) = sv1;
	369	s2.Column(kc) = sv2;
	370	v12re.Column(kc) = real(vv12);
	371	v12im.Column(kc) = imag(vv12);
	372	v12phi.Column(kc) = phase(vv12);
[3647]	373	v12mod.Column(kc) = module(vv12);
	374
	375	// Calcul moyenne dans des bandes en frequence
	376	int deltajf=(jf2-jf1)/nfreq;
	377	if (deltajf<1) deltajf=1;
	378	for(int kf=0; kf<nfreq; kf++) {
	379	sa_size_t jfstart=jf1+kf*deltajf;
	380	sa_size_t jfend=jfstart+deltajf;
	381	if (jfend>jf2) break;
	382	nmnt=0; ms1=ms2=ms12=ms12re=ms12im=ms12phi=ms12mod=0.;
	383	for(sa_size_t jf=jfstart; jf<jfend; jf++) {
	384	ms1 += s1(jf,kc);
	385	ms2 += s2(jf,kc);
	386	ms12re += v12re(jf,kc);
	387	ms12im += v12im(jf,kc);
	388	ms12phi += v12phi(jf,kc);
	389	ms12mod += v12mod(jf,kc);
	390	}
	391	nmnt = (jfend-jfstart);
	392	if (nmnt>0) {
	393	double fnorm = (double)nmnt;
	394	xnt[0] = ((int_8)(sv1.Info()["StartFC"])+(int_8)(sv1.Info()["EndFC"]))*0.5;
	395	xnt[1] = ((int_8)(sv1.Info()["StartTT"])+(int_8)(sv1.Info()["EndTT"]))*0.5;
	396	xnt[2] = kf;
	397	xnt[3] = ms1/fnorm;
	398	xnt[4] = ms2/fnorm;
	399	xnt[5] = ms12/fnorm;
	400	xnt[6] = ms12re/fnorm;
	401	xnt[7] = ms12im/fnorm;
	402	xnt[8] = ms12phi/fnorm;
	403	xnt[9] = ms12mod/fnorm;
	404	nt.Fill(xnt);
	405	}
[3646]	406	}
	407	kc++;
	408
	409	}
	410	if (card==2)
	411	sprintf(fname, "%s/Ch34mtx.ppf",inoutpath.c_str());
	412	else
	413	sprintf(fname, "%s/Ch12mtx.ppf",inoutpath.c_str());
	414
	415	cout << nt;
[3699]	416	cout << "ProcSVFilesVJun09: Opening file " << fname << " for writing" << endl;
[3646]	417	POutPersist po(fname);
	418	string tag1="s1";
	419	string tag2="s2";
	420	string tag12r="v12re";
	421	string tag12i="v12im";
	422	string tag12p="v12phi";
	423	string tagnt="nt12";
	424	if (card==2) {
	425	tag1="s3";
	426	tag2="s4";
	427	tag12r="v34re";
	428	tag12i="v34im";
	429	tag12p="v34phi";
	430	tagnt="nt34";
	431	}
	432	po << PPFNameTag(tag1) << s1;
	433	po << PPFNameTag(tag2) << s2;
	434	po << PPFNameTag(tag12r) << v12re;
	435	po << PPFNameTag(tag12i) << v12im;
	436	po << PPFNameTag(tag12p) << v12phi;
	437	po << PPFNameTag(tagnt) << nt;
[3699]	438	cout << "ProcSVFilesVJun09: Matrices s1, s2, v12re, v12im, v12phi, NTuple nt written to file " << fname << endl;
[3646]	439	return 0;
	440	}
	441
	442

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: Sophya/trunk/AddOn/TAcq/svv2mtx.cc@ 4059

Download in other formats: