source: Sophya/trunk/Cosmo/SimLSS/cmvconcherr.cc@ 3756

Last change on this file since 3756 was 3330, checked in by cmv, 18 years ago

intro ComputeSpectrum avec soustraction de bruit et deconvolution par la fct de transfert du pixel cmv 01/10/2007

File size: 10.5 KB
Line 
1#include "sopnamsp.h"
2#include "machdefs.h"
3#include <iostream>
4#include <stdlib.h>
5#include <stdio.h>
6#include <string.h>
7#include <math.h>
8#include <unistd.h>
9#include "timing.h"
10#include "histerr.h"
11#include "hist2err.h"
12#include "fitshisterr.h"
13
14int ObjectType(string nameh,string nameppf);
15int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp=2,bool do_cov=false);
16int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp=2,int ibinkt=-1,int ibinkz=-1);
17void usage(void);
18void usage(void)
19{
20cout
21<<"cmvconcherr -w wtyp -n name_histoerr file1.ppf file2.ppf ..."<<endl
22<<" name_histoerr: object name"<<endl
23<<" wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl
24<<" -C : compute covariance for 1D spectrum"<<endl
25<<" -Z ibinKt : compute covariance of line Pk(ibinKt,.) along Kz for 2D"<<endl
26<<" -T ibinKz : compute covariance of line Pk(.,ibinKz) along Kt for 2D"<<endl
27<<endl;
28}
29
30
31//---------------------------------------------------------------
32int main(int narg,char *arg[])
33{
34 // "hpkgen" "hpkgenfb" "hpkgenf" "hpkrecb" "hpkrec"
35 // "hpkgen2" "hpkgenfb2" "hpkgenf2" "hpkrecb2" "hpkrec2"
36 string nameh = "";
37 vector<string> ppfname;
38 int Write_Type = 2;
39 bool Do_Cov = false;
40 int ibinKt=-1, ibinKz=-1;
41
42 // --- Decodage arguments
43 char c;
44 while((c = getopt(narg,arg,"hCn:w:Z:T:")) != -1) {
45 switch (c) {
46 case 'n' :
47 nameh = optarg;
48 break;
49 case 'w' :
50 sscanf(optarg,"%d",&Write_Type);
51 break;
52 case 'C' :
53 Do_Cov = true;
54 break;
55 case 'Z' :
56 sscanf(optarg,"%d",&ibinKt);
57 break;
58 case 'T' :
59 sscanf(optarg,"%d",&ibinKz);
60 break;
61 case 'h' :
62 default :
63 usage(); return -1;
64 }
65 }
66 if(nameh.size()<=0 || optind>=narg) {usage(); return -1;}
67 for (int i=optind;i<narg;i++) ppfname.push_back(arg[i]);
68 cout<<"nameh = "<<nameh
69 <<" , number of file to be read = "<<ppfname.size()
70 <<" , write_type="<<Write_Type<<" , cov="<<Do_Cov<<endl;
71
72 // --- Quel type d'objet ?
73 int obtyp = ObjectType(nameh,ppfname[0]);
74
75 // --- lecture et moyenne des HistoErr ou Histo2DErr
76 int nread=0;
77 if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type,Do_Cov);
78 else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type,ibinKt,ibinKz);
79 else {
80 cout<<"Type of object "<<nameh<<" not decoded"<<endl;
81 return -1;
82 }
83 cout<<"Number of file read: "<<nread<<endl;
84 return 0;
85}
86
87//---------------------------------------------------------------
88int ObjectType(string nameh,string nameppf)
89// retourne: 1 si HistoErr
90// 2 si Histo2DErr,
91// 0 si autre objet
92// -1 si il n'y a pas d'objet de nom "nameh"
93{
94 PInPersist pis(nameppf.c_str());
95
96 bool found_tag = pis.GotoNameTag(nameh.c_str());
97 if(!found_tag) return -1;
98
99 PPersist *pps = pis.ReadObject();
100 AnyDataObj *obj = pps->DataObj();
101 cout<<"Object type read from input PPF file : "<<typeid(*obj).name()<< endl;
102
103 HistoErr *herr = dynamic_cast<HistoErr *>(obj);
104 if(herr) return 1;
105
106 Histo2DErr *herr2 = dynamic_cast<Histo2DErr *>(obj);
107 if(herr2) return 2;
108
109 return 0;
110}
111
112//---------------------------------------------------------------
113int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp,bool do_cov)
114{
115 if(ppfname.size()<=0) return 0;
116
117 HistoErr *herrconc = NULL;
118 TVector<r_8> Tsum; TMatrix<r_8> Tsum2;
119 int nherr=0, nread=0, itest=0;
120 double sum=0., sum2=0., nsum=0;
121 for (unsigned int ifile=0;ifile<ppfname.size();ifile++) {
122 PInPersist pis(ppfname[ifile].c_str());
123 HistoErr herr;
124 pis.GetObject(herr,nameh);
125 cout<<ifile<<" : "<<ppfname[ifile]<<" nbin="<<herr.NBins()<<endl;
126
127 if(ifile==0) {
128 herrconc = new HistoErr(herr.XMin(),herr.XMax(),herr.NBins());
129 nherr = herr.NBins();
130 itest = int(herrconc->NBins()/5.+0.5);
131 if(do_cov) {Tsum.ReSize(nherr); Tsum=0.; Tsum2.ReSize(nherr,nherr); Tsum2=0.;}
132 } else if(nherr!=herr.NBins()) {
133 cout<<"BAD NUMBER OF BINS"<<endl;
134 continue;
135 }
136
137 if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
138 <<" nmean="<<herr.NMean()<<endl;
139 nread++;
140 for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
141 sum+=herr(itest); sum2+=herr(itest)*herr(itest); nsum++;
142
143 if(do_cov) {
144 for(int i=0;i<nherr;i++) {
145 Tsum(i) += herr(i);
146 for(int j=0;j<nherr;j++) Tsum2(i,j) += herr(i)*herr(j);
147 }
148 }
149
150 }
151 herrconc->ToVariance();
152 if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
153 cout<<"Test bin "<<itest<<" mean="<<sum<<" sigma^2="<<sum2<<" nsum="<<nsum<<endl;
154 cout<<" "<<" mean="<<(*herrconc)(itest)
155 <<" sigma^2="<<herrconc->Error2(itest)<<endl;
156
157 // --- Covariance
158 if(do_cov && nread>1) {
159 Tsum /= nread;
160 Tsum2 /= nread;
161 for(int i=0;i<nherr;i++)
162 for(int j=0;j<nherr;j++) Tsum2(i,j) -= Tsum(i)*Tsum(j);
163 }
164
165 // --- ecriture
166 if(wrtyp&1) { // ecriture ASCII
167 string asname = "cmvconcherr.data";
168 herrconc->WriteASCII(asname);
169 }
170 if(wrtyp&2 || do_cov) { // ecriture PPF
171 string tagobs = "cmvconcherr.ppf";
172 POutPersist posobs(tagobs);
173 tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
174 if(do_cov) {
175 tagobs = "mean"; posobs.PutObject(Tsum,tagobs);
176 tagobs = "cov"; posobs.PutObject(Tsum2,tagobs);
177 }
178 }
179 if(wrtyp&4) { // ecriture FITS
180 FitsInOutFile fio("!cmvconcherr.fits",FitsInOutFile::Fits_Create);
181 fio << *herrconc;
182 }
183
184 delete herrconc;
185
186 return nread;
187}
188
189//---------------------------------------------------------------
190int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp,int ibinkt,int ibinkz)
191{
192 if(ppfname.size()<=0) return 0;
193
194 Histo2DErr *herrconc = NULL;
195 TVector<r_8> Tsum_kz; TMatrix<r_8> Tsum2_kz;
196 TVector<r_8> Tsum_kt; TMatrix<r_8> Tsum2_kt;
197 int nherrx=0, nherry=0, nread=0, itestx=0, itesty=0;
198 double sum=0., sum2=0., nsum=0;
199 for (unsigned int ifile=0;ifile<ppfname.size();ifile++) {
200 PInPersist pis(ppfname[ifile].c_str());
201 Histo2DErr herr;
202 pis.GetObject(herr,nameh);
203 cout<<ifile<<" : "<<ppfname[ifile]<<" nbin="
204 <<herr.NBinX()<<","<<herr.NBinY()<<endl;
205
206 if(ifile==0) {
207 herrconc = new Histo2DErr(herr.XMin(),herr.XMax(),herr.NBinX()
208 ,herr.YMin(),herr.YMax(),herr.NBinY());
209 nherrx = herr.NBinX(); nherry = herr.NBinY();
210 itestx = int(herrconc->NBinX()/5.+0.5);
211 itesty = int(herrconc->NBinY()/5.+0.5);
212 if(ibinkt>=0) { // matrix de covariance pour une ligne Kz
213 if(ibinkt>=nherrx) ibinkt = nherrx-1;
214 Tsum_kz.ReSize(nherry); Tsum_kz=0.;
215 Tsum2_kz.ReSize(nherry,nherry); Tsum2_kz=0.;
216 }
217 if(ibinkz>=0) { // matrix de covariance pour une ligne Kt
218 if(ibinkz>=nherry) ibinkz = nherry-1;
219 Tsum_kt.ReSize(nherrx); Tsum_kt=0.;
220 Tsum2_kt.ReSize(nherrx,nherrx); Tsum2_kt=0.;
221 }
222 } else if(nherrx!=herr.NBinX() || nherry!=herr.NBinY()) {
223 cout<<"BAD NUMBER OF BINS"<<endl;
224 continue;
225 }
226
227 if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
228 <<" nmean="<<herr.NMean()<<endl;
229 nread++;
230 for(int i=0;i<nherrx;i++)
231 for(int j=0;j<nherry;j++) herrconc->AddBin(i,j,herr(i,j));
232 sum+=herr(itestx,itesty); sum2+=herr(itestx,itesty)*herr(itestx,itesty); nsum++;
233
234 if(ibinkt>=0) {
235 for(int i=0;i<Tsum_kz.Size();i++) {
236 Tsum_kz(i) += herr(ibinkt,i);
237 for(int j=0;j<Tsum_kz.Size();j++) Tsum2_kz(i,j) += herr(ibinkt,i)*herr(ibinkt,j);
238 }
239 }
240
241 if(ibinkz>=0) {
242 for(int i=0;i<Tsum_kt.Size();i++) {
243 Tsum_kt(i) += herr(i,ibinkz);
244 for(int j=0;j<Tsum_kt.Size();j++) Tsum2_kt(i,j) += herr(i,ibinkz)*herr(j,ibinkz);
245 }
246 }
247
248 }
249 herrconc->ToVariance();
250 if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
251 cout<<"Test bin "<<itestx<<","<<itesty
252 <<" mean="<<sum<<" sigma^2="<<sum2<<" nsum="<<nsum<<endl;
253 cout<<" "<<" mean="<<(*herrconc)(itestx,itesty)
254 <<" sigma^2="<<herrconc->Error2(itestx,itesty)<<endl;
255
256 // --- Covariance
257 if(ibinkt>=0 && nread>1) {
258 Tsum_kz /= nread;
259 Tsum2_kz /= nread;
260 for(int i=0;i<Tsum_kz.Size();i++)
261 for(int j=0;j<Tsum_kz.Size();j++) Tsum2_kz(i,j) -= Tsum_kz(i)*Tsum_kz(j);
262 }
263 if(ibinkz>=0 && nread>1) {
264 Tsum_kt /= nread;
265 Tsum2_kt /= nread;
266 for(int i=0;i<Tsum_kt.Size();i++)
267 for(int j=0;j<Tsum_kt.Size();j++) Tsum2_kt(i,j) -= Tsum_kt(i)*Tsum_kt(j);
268 }
269
270 // --- ecriture
271 if(wrtyp&1) { // ecriture ASCII
272 string asname = "cmvconcherr2.data";
273 herrconc->WriteASCII(asname);
274 }
275 if(wrtyp&2) { // ecriture PPF
276 string tagobs = "cmvconcherr2.ppf";
277 POutPersist posobs(tagobs);
278 tagobs = "herrconc2"; posobs.PutObject(*herrconc,tagobs);
279 if(ibinkt>=0) {
280 char str[32];
281 sprintf(str,"kzm%d",ibinkt);
282 posobs.PutObject(Tsum_kz,str);
283 sprintf(str,"kzc%d",ibinkt);
284 posobs.PutObject(Tsum2_kz,str);
285 }
286 if(ibinkz>=0) {
287 char str[32];
288 sprintf(str,"ktm%d",ibinkz);
289 posobs.PutObject(Tsum_kt,str);
290 sprintf(str,"ktc%d",ibinkz);
291 posobs.PutObject(Tsum2_kt,str);
292 }
293 }
294 if(wrtyp&4) { // ecriture FITS
295 FitsInOutFile fio("!cmvconcherr2.fits",FitsInOutFile::Fits_Create);
296 fio << *herrconc;
297 }
298
299 delete herrconc;
300
301 return nread;
302}
303
304/*
305####
306#### HistoErr 1D
307####
308openppf cmvconcherr.ppf
309
310c++exec int i = 700; \
311cout<<i<<" "<<herrconc.BinCenter(i)<<" "<<herrconc(i)<<" "<<herrconc.Error2(i)<<" "<<herrconc.NEntBin(i)<<endl;
312
313zone 2 2
314disp herrconc "hbincont err"
315disp herrconc "hbinerr"
316disp herrconc "hbinent"
317
318zone
319n/plot herrconc.val%log10(x) x>0 ! "connectpoints"
320n/plot herrconc.val-sqrt(err2)%log10(x) x>0&&err2>0 ! "connectpoints same red"
321n/plot herrconc.val+sqrt(err2)%log10(x) x>0&&err2>0 ! "connectpoints same red"
322
323n/plot herrconc.sqrt(err2)%log10(x) x>0&&err2>0 ! "connectpoints same red"
324n/plot herrconc.sqrt(err2)/val%log10(x) x>0&&err2>0&&val>0 ! "connectpoints"
325
326disp mean
327imag cov
328
329del cor
330c++exec \
331TMatrix<r_8> cor(cov,false); cor = 0.; \
332for(int i=0;i<cor.NRows();i++) { \
333 for(int j=0;j<cor.NCols();j++) { \
334 double v = cov(i,i)*cov(j,j); \
335 if(v<=0.) continue; \
336 cor(i,j) = cov(i,j)/sqrt(v); \
337} } \
338KeepObj(cor); \
339cout<<"Matrice cor cree "<<endl;
340
341imag cor
342n/plot cor.val%c r==500 ! "connectpoints"
343
344####
345#### Histo2DErr 2D
346####
347openppf cmvconcherr2.ppf
348
349zone 2 2
350imag herrconc2 "hbincont"
351imag herrconc2 "hbinerr"
352imag herrconc2 "hbinent"
353
354zone
355
356set mean kzm140
357set cov kzc140
358
359set mean ktm32
360set cov ktc32
361
362n/plot $mean.val%n ! ! "connectpoints"
363disp $cov
364n/plot $cov.val%c r==32 ! "connectpoints"
365
366set cor ${cov}cor
367del $cor
368c++exec \
369TMatrix<r_8> $cor($cov,false); $cor = 0.; \
370for(int i=0;i<$cor.NRows();i++) { \
371 for(int j=0;j<$cor.NCols();j++) { \
372 double v = $cov(i,i)*$cov(j,j); \
373 if(v<=0.) continue; \
374 $cor(i,j) = $cov(i,j)/sqrt(v); \
375} } \
376KeepObj($cor); \
377cout<<"Matrice cor cree "<<endl;
378
379disp $cor
380n/plot $cor.val%c r==32 ! "connectpoints"
381*/
Note: See TracBrowser for help on using the repository browser.