source: Sophya/trunk/ArchTOIPipe/ProcWSophya/toi2map.cc@ 2227

Last change on this file since 2227 was 2160, checked in by cmv, 23 years ago

BUG constructeur localmap par copie En shared cvstatcvstat cmv 5/8/2002
File size: 7.9 KB
Line 
1// ArchTOIPipe (C) CEA/DAPNIA/SPP IN2P3/LAL
2// Eric Aubourg
3// Christophe Magneville
4// Reza Ansari
5// $Id: toi2map.cc,v 1.30 2002-08-05 16:25:13 cmv Exp $
6
7#include "machdefs.h"
8#include "toimanager.h"
9#include "pexceptions.h"
10#include "ctimer.h"
11#include "toi2map.h"
12// La valeur "Pi" doit etre celle de smathconst.h a cause du test sur theta
13#include "smathconst.h"
14
15////////////////////////////////////////////////////////////////////////
16TOI2Map::TOI2Map(PixelMap<r_8>* map,PixelMap<r_8>* wmap)
17 : mMap(map), mWMap(wmap), mWMapInternal(false), totnscount(0), usednscount(0)
18{
19 SetEquinox();
20 SetCoorIn();
21 SetCoorMap();
22 SetCalibrationFactor();
23 SetTestFlag();
24 SetTestMin();
25 SetTestMax();
26
27 if(mMap->NbPixels()<1) {
28 cout<<"TOI2Map::TOI2Map() Bad number of pixels in sphere mMap "
29 <<mMap->NbPixels()<<endl;
30 throw ParmError("TOI2Map::TOI2Map() - Bad number of pixels in sphere");
31 }
32 mMap->SetPixels(0.);
33 int nlat = mMap->SizeIndex();
34 string typmap = mMap->TypeOfMap();
35
36 if(mWMap==NULL) {
37 // We would need a cloning function in maps. $$TBD$$ $CHECK$
38 if (typmap == "LOCAL") {
39 mWMap = new LocalMap<r_8>(*(LocalMap<r_8>*)mMap,false);
40 } else if (typmap == "RING") {
41 mWMap = new SphereHEALPix<r_8>(nlat);
42 } else {
43 cout << "TOI2Map::TOI2Map() cannot handle map of type " << typmap << endl;
44 throw ParmError("TOI2Map::TOI2Map() - bad type of map");
45 }
46 mWMapInternal = true;
47 } else {
48 mWMapInternal = false;
49 if(nlat != mWMap->SizeIndex()) {
50 cout<<"TOI2Map::TOI2Map() Bad size for sphere mWMap, does not "
51 <<"correspond to mMap : "<<mMap->SizeIndex()<<", "<<mWMap->SizeIndex()<<endl;
52 if (typmap == "LOCAL") {
53 throw ParmError("TOI2Map::TOI2Map() - Different sizes for map and wmap");
54 } else if(typmap == "RING") {
55 ((SphereHEALPix<r_8> *)mWMap)->Resize(nlat);
56 cout<<"Resize have been done..."<<endl;
57 } else {
58 cout << "TOI2Map::TOI2Map() cannot handle map of type " << typmap << endl;
59 throw ParmError("TOI2Map::TOI2Map() - bad type of map");
60 }
61 }
62 }
63 if(mWMap->NbPixels()<1) {
64 cout<<"TOI2Map::TOI2Map() Bad number of pixels in sphere mWMap "
65 <<mWMap->NbPixels()<<endl;
66 throw ParmError("TOI2Map::TOI2Map() - Bad number of pixels in sphere");
67 }
68 mWMap->SetPixels(0);
69}
70
71TOI2Map::~TOI2Map()
72{
73 if(mWMap && mWMapInternal) delete mWMap;
74}
75
76////////////////////////////////////////////////////////////////////////
77void TOI2Map::Print(::ostream & os)
78{
79 os<<"TOI2Map::Print -- Map type " << mMap->TypeOfMap() << " SizeIndex = "<<mMap->SizeIndex()<<endl
80 <<" - Equinoxe="<<mActualYear<<endl
81 <<" - TypCoorIn: "<<mTypCoorIn<<" = "<<DecodeTypAstro(mTypCoorIn)<<endl
82 <<" - TypCoorMap: "<<mTypCoorMap<<" = "<<DecodeTypAstro(mTypCoorMap)<<endl
83 <<" - Tests: Flag("<<mTFlag<<") bad="<<mBadFlag
84 <<" / Value Min("<<mTMin<<")="<<mValMin
85 <<" , Max("<<mTMax<<")="<<mValMax<<endl;
86 os<<" TotalProcessedSamples="<<ProcessedSampleCount()
87 <<" UsedSampleCount="<<UsedSampleCount()<<endl;
88}
89
90////////////////////////////////////////////////////////////////////////
91void TOI2Map::init() {
92 cout << "TOI2Map::init" << endl;
93 declareInput("Coord1In"); // input index 0
94 declareInput("Coord2In"); // input index 1
95 declareInput("BoloIn"); // input index 2
96}
97
98////////////////////////////////////////////////////////////////////////
99// define SANS_BUFFER
100void TOI2Map::run()
101{
102long snb = getMinIn();
103long sne = getMaxIn();
104
105if(snb>sne) {
106 cout<<"TOI2Map::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
107 throw ParmError("TOI2Map::run() - Bad sample interval");
108}
109if(!checkInputTOIIndex(0) || !checkInputTOIIndex(1) || !checkInputTOIIndex(2)) {
110 cout<<"TOI2Map::run() - Input TOI (Coord1In or Coord2In or BoloIn) not connected! "<<endl;
111 throw ParmError("TOI2Map::run() Output TOI (Coord1In or Coord2In or BoloIn) not connected!");
112}
113if( !(mTypCoorIn&TypCoordEq || mTypCoorIn&TypCoordGal) ) {
114 cout<<"TOI2Map::run() - Input Coordinates not Eq or Gal! "<<endl;
115 throw ParmError("TOI2Map::run() - Input Coordinates not Eq or Gal!");
116}
117if( !(mTypCoorMap&TypCoordEq || mTypCoorMap&TypCoordGal) ) {
118 cout<<"TOI2Map::run() - Output Coordinates not Eq or Gal! "<<endl;
119 throw ParmError("TOI2Map::run() - Output Coordinates not Eq or Gal!");
120}
121
122//---------------------------------------------------------
123#define NFILL 25
124try {
125
126int ii;
127uint_4 mNSnFill=0, mNpixFill=0, NFill[NFILL], BadCoorRange=0, nbadsn=0;
128for(ii=0;ii<NFILL;ii++) NFill[ii]=0;
129double mjd = MJDfrYear(mActualYear);
130
131cout<<"TOI2Map::run() from "<<snb<<" to "<<sne;
132#ifndef SANS_BUFFER
133int bufsz = 100;
134uint_8* bfgbolo = new uint_8[bufsz];
135double* bbolo = new double[bufsz];
136double* bc1 = new double[bufsz];
137double* bc2 = new double[bufsz];
138int i0 = -1;
139cout<<" (getData() bufferized)"<<endl;
140#else
141cout<<" (getData() not bufferized)"<<endl;
142#endif
143
144// Remplissage des spheres
145for(int s=snb;s<=sne;s++) {
146 uint_8 fgbolo = 0;
147 double bolo,coord1,coord2;
148 // Equatoriales / Galactiques
149 // coord1,2 = alpha,delta / gLon,gLat
150
151#ifndef SANS_BUFFER
152 if(i0<0 || s<i0 || s>=i0+bufsz) {
153 i0 = s;
154 int nget = (sne-s+1<bufsz)? nget=sne-s+1: bufsz;
155 getData(0,i0,nget,bc1);
156 getData(1,i0,nget,bc2);
157 getData(2,i0,nget,bbolo,bfgbolo);
158 totnscount += nget;
159 }
160 bolo = bbolo[s-i0]; fgbolo = bfgbolo[s-i0];
161 coord1 = bc1[s-i0]; coord2 = bc2[s-i0];
162#else
163 getData(2,s,bolo,fgbolo);
164 coord1 = getData(0,s);
165 coord2 = getData(1,s);
166 totnscount++;
167#endif
168
169 // Comme il n'y a pas de toi en sortie, il faut dire
170 // aux processeur/toi que l'on a plus besoin des donnees.
171 if (s%100 == 0) wontNeedBefore(s-1);
172 //if(fgbolo>0) cout<<fgbolo<<" , "<<(fgbolo&mBadFlag)<<endl;
173 if(mTFlag && (fgbolo&mBadFlag)) {nbadsn++; continue;}
174 if(mTMin && (bolo<mValMin)) {nbadsn++; continue;}
175 if(mTMax && (bolo>mValMax)) {nbadsn++; continue;}
176 usednscount++;
177
178 // sphere phi entre [0,2*Pi] en radians
179 // sphere theta entre [0,Pi] en radians
180 double phi=-1.;
181 CoordConvertToStd(mTypCoorIn,&coord1,&coord2);
182
183 if(mTypCoorIn&TypCoordEq && mTypCoorMap&TypCoordGal) { // Eq -> Gal
184 EqtoGal(mjd,coord1,coord2,&coord1,&coord2);
185 phi = coord1 * Pi/180.;
186 } else if(mTypCoorIn&TypCoordGal && mTypCoorMap&TypCoordEq) { // Gal -> Eq
187 GaltoEq(mjd,coord1,coord2,&coord1,&coord2);
188 phi = coord1 * Pi/12.;
189 } else if(mTypCoorMap&TypCoordGal) { // Gal -> Gal
190 phi = coord1 * Pi/180.;
191 } else if(mTypCoorMap&TypCoordEq) { // Eq -> Eq
192 phi = coord1 * Pi/12.;
193 }
194 ToCoLat(&coord2,TypUniteD);
195 double theta = coord2 * Pi/180.;
196
197 if(phi<0. | phi>=2*Pi || theta<0. || theta>Pi)
198 {BadCoorRange++; continue;}
199
200 int_4 ipix = mMap->PixIndexSph(theta,phi);
201 if ((ipix < 0) || (ipix >= mMap->NbPixels()) ) continue;
202 (*mMap)(ipix) += bolo*mCalibFactor;
203 ((*mWMap)(ipix)) += 1;
204 mNSnFill++;
205}
206
207 cout<<"TOI2Map::run(): End of SN loop - TotalProcessedSamples="
208 <<ProcessedSampleCount()<<" UsedSampleCount="<<UsedSampleCount()<<endl;
209
210// Remplissage des spheres
211 for(int_4 i=0;i<mMap->NbPixels();i++) {
212 r_8 wf = (*mWMap)(i);
213 if(wf>0.) {mNpixFill++; (*mMap)(i) /= wf;}
214 int_4 nf = int_4(wf);
215 if(nf>=NFILL) nf=NFILL-1; NFill[nf]++;
216 }
217
218 cout<<"TOI2Map::run(): mNpixTot="<<mMap->NbPixels()
219 <<" mNpixFill="<<mNpixFill
220 <<" --> FracSky="<<mNpixFill*100./(double)mMap->NbPixels()<<"%"
221 <<" NFill["<<NFILL<<"] ="<<endl;
222 for(ii=0;ii<NFILL;ii++) {cout<<NFill[ii]<<" "; if(ii%10==9) cout<<endl;}
223 cout<<endl;
224 cout<<" mNSnFill="<<mNSnFill<<" NBadSn="<<nbadsn
225 <<" BadCoorRange="<<BadCoorRange<<endl;
226
227#ifndef SANS_BUFFER
228delete [] bbolo; delete [] bfgbolo;
229delete [] bc1; delete [] bc2;
230#endif
231
232//---------------------------------------------------------
233} catch (PException & exc) {
234 cout<<"TOI2Map: Catched Exception "<<(string)typeid(exc).name()
235 <<"\n .... Msg= "<<exc.Msg()<<endl;
236}
237
238return;
239}
Note: See TracBrowser for help on using the repository browser.