/[PAMELA software]/quicklook/dataToXML/OrbitalRate.cpp
ViewVC logotype

Annotation of /quicklook/dataToXML/OrbitalRate.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (hide annotations) (download)
Wed May 31 07:11:04 2006 UTC (18 years, 6 months ago) by kusanagi
Branch: MAIN
CVS Tags: dataToXML1_02/00, dataToXML1_01/00
Changes since 1.3: +7 -3 lines
Bug Fix.
- Now check if the mcmdItem is not NULL.

1 kusanagi 1.1 /**
2     * Rate
3     * author Nagni
4     * version 1.0 - 27 April 2006
5     *
6     * Description: .
7     *
8     * Parameters:
9     * base - the path where to find the PAMELA unpacked root file.
10     * outDirectory - the path where to put the output file.
11     * orbPath - the path where to find an orb format for the output.
12     *
13     * version 1.0 - 27 April 2006
14     * First implementation
15     */
16    
17     #include <mcmd/McmdEvent.h>
18     #include <mcmd/McmdRecord.h>
19     #include <EventHeader.h>
20     #include <PscuHeader.h>
21     #include <TTree.h>
22     #include "cTle.h"
23     #include "cEci.h"
24     #include "cOrbit.h"
25     #include "TH2F.h"
26     #include "TFrame.h"
27     #include "TGraph.h"
28     #include "TCanvas.h"
29     #include "TASImage.h"
30     #include "TMarker.h"
31    
32     #include "TString.h"
33     #include "TObjString.h"
34     #include "TStyle.h"
35     #include "TPaletteAxis.h"
36     #include "TROOT.h"
37 kusanagi 1.2 #include <sys/stat.h>
38 kusanagi 1.3 #include <fstream>
39 kusanagi 1.1
40     #define TRUE 1
41     #define FALSE 0
42     /**
43     *
44     * @param base
45     * @param outDirectory
46     * @param xslPath
47     */
48 kusanagi 1.2 using namespace std;
49 kusanagi 1.1
50     void InitStyle() {
51     gROOT->SetStyle("Plain");
52    
53     TStyle *myStyle[2], *tempo;
54     myStyle[0]=new TStyle("StyleWhite", "white");
55     myStyle[1]=new TStyle("StyleBlack", "black");
56    
57     tempo=gStyle;
58     Int_t linecol, bkgndcol, histcol;
59    
60     for(Int_t style=0; style<2; style++) {
61    
62     linecol=kWhite*style+kBlack*(1-style);
63     bkgndcol=kBlack*style+kWhite*(1-style);
64     histcol=kYellow*style+kBlack*(1-style); // was 95
65    
66     myStyle[style]->Copy(*tempo);
67    
68     myStyle[style]->SetCanvasBorderMode(0);
69     myStyle[style]->SetCanvasBorderSize(1);
70     myStyle[style]->SetFrameBorderSize(1);
71     myStyle[style]->SetFrameBorderMode(0);
72     myStyle[style]->SetPadBorderSize(1);
73     myStyle[style]->SetStatBorderSize(1);
74     myStyle[style]->SetTitleBorderSize(1);
75     myStyle[style]->SetPadBorderMode(0);
76     myStyle[style]->SetPalette(1,0);
77     myStyle[style]->SetPaperSize(20,27);
78     myStyle[style]->SetFuncColor(kRed);
79     myStyle[style]->SetFuncWidth(1);
80     myStyle[style]->SetLineScalePS(1);
81     myStyle[style]->SetCanvasColor(bkgndcol);
82     myStyle[style]->SetAxisColor(linecol,"XYZ");
83     myStyle[style]->SetFrameFillColor(bkgndcol);
84     myStyle[style]->SetFrameLineColor(linecol);
85     myStyle[style]->SetLabelColor(linecol,"XYZ");
86     myStyle[style]->SetPadColor(bkgndcol);
87     myStyle[style]->SetStatColor(bkgndcol);
88     myStyle[style]->SetStatTextColor(linecol);
89     myStyle[style]->SetTitleColor(linecol,"XYZ");
90     myStyle[style]->SetTitleFillColor(bkgndcol);
91     myStyle[style]->SetTitleTextColor(linecol);
92     myStyle[style]->SetLineColor(linecol);
93     myStyle[style]->SetMarkerColor(histcol);
94     myStyle[style]->SetTextColor(linecol);
95    
96     myStyle[style]->SetGridColor((style)?13:kBlack);
97     myStyle[style]->SetHistFillStyle(1001*(1-style));
98     myStyle[style]->SetHistLineColor(histcol);
99     myStyle[style]->SetHistFillColor((style)?bkgndcol:kYellow);
100     }
101    
102     myStyle[1]->cd();
103    
104     gROOT->ForceStyle();
105    
106     }
107    
108 kusanagi 1.3 void Rate(TString base, TString outDirectory = "", TString format = "jpg", TString mapFile = "", TString tleFile = ""){
109 kusanagi 1.1 TTree *tr = 0;
110     pamela::McmdEvent *mcmdev = 0;
111     pamela::McmdRecord *mcmdrc = 0;
112     TArrayC *mcmddata = 0;
113     ULong64_t nevents = 0;
114     ULong64_t timesync = 0;
115     pamela::EventHeader *eh = 0;
116     pamela::PscuHeader *ph = 0;
117     stringstream oss;
118     double offsetTime = 0;
119     double absTime;
120     UInt_t i = 0;
121     UInt_t j = 0;
122 kusanagi 1.2 struct stat buf;
123 kusanagi 1.1
124 kusanagi 1.2 i = stat (mapFile.Data(), &buf );
125     // If the file does not exists
126     if (i != 0){
127     printf("The %s file does not exists.", mapFile.Data());
128     exit(0);
129     }
130 kusanagi 1.1
131     TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
132     filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();
133     // Test SGP4
134 kusanagi 1.3 string str1 = "SGP4 Test";
135 kusanagi 1.1 string str2 = "1 25544U 98067A 06117.32388940 .00009459 00000-0 64427-4 0 8131";
136     string str3 = "2 25544 51.6388 89.2928 0009043 155.3293 354.6512 15.75673697425143";
137 kusanagi 1.3 if (tleFile != ""){
138     fstream fileTle(tleFile.Data(),ios::in);
139     if (fileTle.is_open()) {
140     getline (fileTle,str1);
141     getline (fileTle,str2);
142     getline (fileTle,str3);
143     }
144     fileTle.close();
145     }
146 kusanagi 1.1
147     cTle tle1(str1, str2, str3);
148     cOrbit orbit(tle1);
149     cEci eci;
150     cCoordGeo coo;
151     TH2F *rate = new TH2F("rate", base, 360, -180, 180, 161, -80.5, 80.5);
152     TFile *rootFile = new TFile(base);
153     if (rootFile->IsZombie()) {
154     printf("The %s file does not exist", base.Data());
155     exit(0);
156     }
157    
158    
159     /*
160     * process Mcmd
161     */
162     long int recEntries;
163     tr = (TTree*)rootFile->Get("Mcmd");
164     tr->SetBranchAddress("Mcmd", &mcmdev);
165     nevents = tr->GetEntries();
166    
167     bool timeFound = FALSE;
168 kusanagi 1.4 while (i < nevents) {
169 kusanagi 1.1 tr->GetEntry(i);
170     recEntries = mcmdev->Records->GetEntries();
171 kusanagi 1.4 while (j < recEntries){
172 kusanagi 1.1 mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);
173 kusanagi 1.4 if (not (mcmdrc == NULL) && mcmdrc->ID1 == 0xE0){
174 kusanagi 1.1 mcmddata = mcmdrc->McmdData;
175     timesync = (((ULong64_t)mcmddata->At(0)<<24)&0xFF000000) +
176     (((ULong64_t)mcmddata->At(1)<<16)&0x00FF0000) +
177     (((ULong64_t)mcmddata->At(2)<<8)&0x0000FF00) +
178     (((ULong64_t)mcmddata->At(3))&0x000000FF);
179     offsetTime = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
180     timeFound = TRUE;
181     }
182     j++;
183     }
184     i++;
185     }
186    
187 kusanagi 1.4 if (!timeFound) {
188     printf("No timesync info have been found in the file %s", base.Data());
189     exit(0);
190     }
191 kusanagi 1.1 tr = (TTree*)rootFile->Get("Physics");
192     TBranch *headBr = tr->GetBranch("Header");
193     tr->SetBranchAddress("Header", &eh);
194     nevents = tr->GetEntries();
195     //Fill variables from root-ple
196     for (i = 0; i < nevents; i++){
197     tr->GetEntry(i);
198     ph = eh->GetPscuHeader();
199     absTime = ((ph->GetOrbitalTime()*(1./1000.)) + absTime)/60;
200     orbit.getPosition(absTime, &eci);
201     coo = eci.toGeo();
202     rate->Fill(rad2deg(coo.m_Lon), rad2deg(coo.m_Lat));
203     /*
204     printf(" %16.8f %16.8f %16.8f\n",
205     rad2deg(coo.m_Lat),
206     rad2deg(coo.m_Lon),
207     coo.m_Alt);
208     */
209     }
210     double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000;
211    
212     int ptcnt=0, color=0;
213     TMarker *tma=NULL;
214     TLine *tli=NULL;
215     double step=0;
216    
217 kusanagi 1.2
218     TImage *tImage=TImage::Open(mapFile);
219 kusanagi 1.1 int width=(int)(tImage->GetWidth()*0.80);
220     int height=(int)(tImage->GetHeight()*0.80);
221     InitStyle();
222     TCanvas *c1 = new TCanvas("c1","rate/orbit",-width, height); // - : removes the menu bars
223    
224     TH1F *hframe=NULL;
225     hframe=gPad->DrawFrame(-180,-90,180,90);
226     oss.str("");
227     oss << filename << " - Event Rate (Hz)";
228     hframe->SetTitle(oss.str().c_str());
229     hframe->SetXTitle("Longitude (deg)");
230     hframe->SetYTitle("Latitude (deg)");
231     c1->cd();
232    
233    
234    
235     TPad *p2 = new TPad("p2", "p2", 0.10, 0.04, 0.983, 1);
236     p2->Draw();
237     p2->cd();
238     TPaletteAxis *tpa=NULL;
239     TH2F *forpal=new TH2F("forpal","",2,0,2,2,0,2);
240     forpal->SetAxisColor(kBlack); //Delete the stat box
241     forpal->SetStats(0); //Delete the stat box
242     forpal->SetMinimum(0.1);
243 kusanagi 1.3 forpal->SetMaximum(15);
244 kusanagi 1.1 forpal->SetBinContent(5,1); // just to initialize the histo
245     forpal->SetContour(50);
246     TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette"));
247     step=forpal->GetContourLevel(1)-forpal->GetContourLevel(0);
248     tpa=new TPaletteAxis(184,-90,195,90,forpal);
249     tpa->SetLabelColor(kWhite);
250     forpal->Draw("colz");
251    
252     c1->cd();
253     TPad *p1 = new TPad("p1", "p1", 0.10,0.10,0.90,0.92);
254     p1->Draw();
255     p1->cd();
256     tImage->Draw("X");
257     c1->cd();
258     gPad->RedrawAxis();
259     c1->Update();
260    
261     c1->cd();
262     ptcnt=0;
263     tma=new TMarker();
264     tma->SetMarkerStyle(4);
265     tli=new TLine();
266     tli->SetLineColor(kMagenta);
267     Stat_t freq;
268     for (Int_t kk = 0; kk < 360; kk++){
269     for (Int_t jj = 0; jj < 161; jj++){
270     freq = rate->GetBinContent(kk, jj);
271     if (freq>0) {
272     posx=(kk - 180); posy= jj - 80.5;
273     // color: palette colors from 51 to 100 ie 50 levels
274     color=51+(int) ((log10((rate==0)?0.1:freq)-log10(0.1))/step); // step defined by palette
275     if (color>100) color=100; // just in case if max rate is not max...
276     tma->SetMarkerColor(color);
277     if (!(posx<0 && oldposx>0) && oldposx!=-1000 && oldposy!=-1000) tli->DrawLine(oldposx,oldposy,posx,posy);
278     tma->DrawMarker(posx,posy);
279     oldposx=posx;
280     oldposy=posy;
281     }
282     }
283     }
284     oss.str("");
285     if (outDirectory == "") {
286     oss << filename.Data() << "_OrbitRate." << format.Data();
287     } else {
288     oss << outDirectory.Data() << filename.Data() << "_OrbitRate." << format.Data();
289     }
290     c1->SaveAs(oss.str().c_str());
291     rootFile->Close();
292     }
293    
294     int main(int argc, char* argv[]){
295 kusanagi 1.2 TString outDir = "./";
296     TString mapFile = "";
297 kusanagi 1.3 TString tleFile = "";
298 kusanagi 1.2 TString format = "jpg";
299 kusanagi 1.1
300     if (argc < 2){
301 kusanagi 1.2 printf("You have to insert at least the file to analyze and the mapFile \n");
302 kusanagi 1.1 printf("Try '--help' for more information. \n");
303     exit(1);
304     }
305    
306     if (!strcmp(argv[1], "--help")){
307 kusanagi 1.2 printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");
308     printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");
309 kusanagi 1.1 printf( "\t --help Print this help and exit \n");
310 kusanagi 1.3 printf( "\t -tleFile[path] Path where to find the tle infos [default dummyOrbit] \n");
311     printf( "\t\tThe tle file have to satisfy a 3-line structure like (this is the included dummyOrbit)\n ");
312     printf( "\t\t\tGP4 Test\n");
313     printf( "\t\t\t1 25544U 98067A 06117.32388940 .00009459 00000-0 64427-4 0 8131\n");
314     printf( "\t\t\t2 25544 51.6388 89.2928 0009043 155.3293 354.6512 15.75673697425143\n");
315 kusanagi 1.1 printf( "\t -outDir[path] Path where to put the output [default ~/tmp] \n");
316     printf( "\t -format[jpg|gif|ps] Format for output files [default 'jpg'] \n");
317     exit(1);
318     }
319    
320     for (int i = 2; i < argc; i++){
321     if (!strcmp(argv[i], "-outDir")){
322     if (++i >= argc){
323     printf( "-outDir needs arguments. \n");
324     printf( "Try '--help' for more information. \n");
325     exit(1);
326     } else {
327     outDir = argv[i];
328     continue;
329     }
330     }
331    
332 kusanagi 1.3 if (!strcmp(argv[i], "-tle")){
333     if (++i >= argc){
334     printf( "-tle needs arguments. \n");
335     printf( "Try '--help' for more information. \n");
336     exit(1);
337     } else {
338     tleFile = argv[i];
339     continue;
340     }
341     }
342    
343 kusanagi 1.2 if (!strcmp(argv[i], "-map")){
344     if (++i >= argc){
345     printf( "-map needs arguments. \n");
346     printf( "Try '--help' for more information. \n");
347     exit(1);
348     } else {
349     mapFile = argv[i];
350     continue;
351     }
352     }
353    
354     if (!strcmp(argv[i], "-format")) {
355 kusanagi 1.1 if (++i >= argc){
356     printf( "-format needs arguments. \n");
357     printf( "Try '--help' for more information. \n");
358     exit(1);
359     } else {
360     format = argv[i];
361     continue;
362     }
363     }
364     }
365 kusanagi 1.2 if (mapFile != ""){
366     Rate(argv[1], outDir, format, mapFile);
367     } else {
368     printf("You have to insert at least the file to analyze and the mapFile \n");
369     printf("Try '--help' for more information. \n");
370     }
371     }
372    
373 kusanagi 1.1
374    
375    

  ViewVC Help
Powered by ViewVC 1.1.23