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

Annotation of /quicklook/dataToXML/OrbitalRate.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (hide annotations) (download)
Fri Dec 8 09:52:00 2006 UTC (17 years, 11 months ago) by pam-rm2
Branch: MAIN
CVS Tags: HEAD
Changes since 1.6: +0 -0 lines
FILE REMOVED
Removed OrbitalRate from this module.

Nico

1 kusanagi 1.1 /**
2     * Rate
3     * author Nagni
4     * version 1.0 - 27 April 2006
5     *
6     * Description: .
7     *
8     *
9 kusanagi 1.6 * 27 April 2006
10 kusanagi 1.1 * First implementation
11 kusanagi 1.6 *
12     * 6 July 2006
13     * Released version adding as an "emergency" Trigger and Anticounter plot. A future release should
14     * foreseen a general class called OrbitalRate, to be put into the "Utility" for plotting whatever
15     * data related to the orbit. Under a preliminary analysis this Class will contain
16     *
17     * OrbitalRate(String TLEFile) - file containing one or more TLE
18     *
19     * void populateHisto(double data, double timeFromTLE) - fill the histogram
20     * TCanvas getCanvas() - return the populated TCanvas
21     *
22 kusanagi 1.1 */
23 kusanagi 1.6 #include <physics/anticounter/AnticounterEvent.h>
24     #include <physics/trigger/TriggerEvent.h>
25 kusanagi 1.1 #include <mcmd/McmdEvent.h>
26     #include <mcmd/McmdRecord.h>
27     #include <EventHeader.h>
28     #include <PscuHeader.h>
29     #include <TTree.h>
30     #include "cTle.h"
31     #include "cEci.h"
32     #include "cOrbit.h"
33 kusanagi 1.5 #include "cJulian.h"
34 kusanagi 1.1 #include "TH2F.h"
35     #include "TFrame.h"
36     #include "TGraph.h"
37     #include "TCanvas.h"
38     #include "TASImage.h"
39     #include "TMarker.h"
40 kusanagi 1.5 #include <TDatime.h>
41 kusanagi 1.1
42     #include "TString.h"
43     #include "TObjString.h"
44     #include "TStyle.h"
45     #include "TPaletteAxis.h"
46     #include "TROOT.h"
47 kusanagi 1.2 #include <sys/stat.h>
48 kusanagi 1.3 #include <fstream>
49 kusanagi 1.1
50     #define TRUE 1
51     #define FALSE 0
52     /**
53     *
54     * @param base
55     * @param outDirectory
56     * @param xslPath
57     */
58 kusanagi 1.2 using namespace std;
59 kusanagi 1.1
60     void InitStyle() {
61     gROOT->SetStyle("Plain");
62    
63     TStyle *myStyle[2], *tempo;
64     myStyle[0]=new TStyle("StyleWhite", "white");
65     myStyle[1]=new TStyle("StyleBlack", "black");
66    
67     tempo=gStyle;
68     Int_t linecol, bkgndcol, histcol;
69    
70     for(Int_t style=0; style<2; style++) {
71    
72     linecol=kWhite*style+kBlack*(1-style);
73     bkgndcol=kBlack*style+kWhite*(1-style);
74     histcol=kYellow*style+kBlack*(1-style); // was 95
75    
76     myStyle[style]->Copy(*tempo);
77    
78     myStyle[style]->SetCanvasBorderMode(0);
79     myStyle[style]->SetCanvasBorderSize(1);
80     myStyle[style]->SetFrameBorderSize(1);
81     myStyle[style]->SetFrameBorderMode(0);
82     myStyle[style]->SetPadBorderSize(1);
83     myStyle[style]->SetStatBorderSize(1);
84     myStyle[style]->SetTitleBorderSize(1);
85     myStyle[style]->SetPadBorderMode(0);
86     myStyle[style]->SetPalette(1,0);
87     myStyle[style]->SetPaperSize(20,27);
88     myStyle[style]->SetFuncColor(kRed);
89     myStyle[style]->SetFuncWidth(1);
90     myStyle[style]->SetLineScalePS(1);
91     myStyle[style]->SetCanvasColor(bkgndcol);
92     myStyle[style]->SetAxisColor(linecol,"XYZ");
93     myStyle[style]->SetFrameFillColor(bkgndcol);
94     myStyle[style]->SetFrameLineColor(linecol);
95     myStyle[style]->SetLabelColor(linecol,"XYZ");
96     myStyle[style]->SetPadColor(bkgndcol);
97     myStyle[style]->SetStatColor(bkgndcol);
98     myStyle[style]->SetStatTextColor(linecol);
99     myStyle[style]->SetTitleColor(linecol,"XYZ");
100     myStyle[style]->SetTitleFillColor(bkgndcol);
101     myStyle[style]->SetTitleTextColor(linecol);
102     myStyle[style]->SetLineColor(linecol);
103     myStyle[style]->SetMarkerColor(histcol);
104     myStyle[style]->SetTextColor(linecol);
105    
106     myStyle[style]->SetGridColor((style)?13:kBlack);
107     myStyle[style]->SetHistFillStyle(1001*(1-style));
108     myStyle[style]->SetHistLineColor(histcol);
109     myStyle[style]->SetHistFillColor((style)?bkgndcol:kYellow);
110     }
111    
112     myStyle[1]->cd();
113    
114     gROOT->ForceStyle();
115    
116     }
117    
118 kusanagi 1.6 void Rate(TString base, TString outDirectory = "", TString format = "png", TString mapFile = "", TString tleFile = "", int offDate = 20060614, int offTime = 210000){
119 kusanagi 1.1 TTree *tr = 0;
120     pamela::McmdEvent *mcmdev = 0;
121     pamela::McmdRecord *mcmdrc = 0;
122 kusanagi 1.5 pamela::EventHeader *eh = 0;
123     pamela::PscuHeader *ph = 0;
124     TArrayC *mcmddata;
125 kusanagi 1.1 ULong64_t nevents = 0;
126 kusanagi 1.5 double timesync = 0;
127 kusanagi 1.1 stringstream oss;
128     double offsetTime = 0;
129 kusanagi 1.5 double absTime = 0;
130 kusanagi 1.1 UInt_t i = 0;
131     UInt_t j = 0;
132 kusanagi 1.2 struct stat buf;
133 kusanagi 1.1
134 kusanagi 1.2 i = stat (mapFile.Data(), &buf );
135     // If the file does not exists
136     if (i != 0){
137     printf("The %s file does not exists.", mapFile.Data());
138     exit(0);
139     }
140 kusanagi 1.1
141     TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
142     filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();
143 kusanagi 1.5
144     string str1 = "RESURS-DK 1";
145 kusanagi 1.6 string str2 = "1 29228U 06021A 06170.19643714 .00009962 00000-0 21000-3 0 196";
146     string str3 = "2 29228 069.9363 054.7893 0167576 127.4359 017.0674 15.31839265 604";
147 kusanagi 1.3 if (tleFile != ""){
148     fstream fileTle(tleFile.Data(),ios::in);
149     if (fileTle.is_open()) {
150     getline (fileTle,str1);
151     getline (fileTle,str2);
152     getline (fileTle,str3);
153     }
154     fileTle.close();
155     }
156 kusanagi 1.1
157     cTle tle1(str1, str2, str3);
158     cOrbit orbit(tle1);
159     cEci eci;
160     cCoordGeo coo;
161     TH2F *rate = new TH2F("rate", base, 360, -180, 180, 161, -80.5, 80.5);
162     TFile *rootFile = new TFile(base);
163     if (rootFile->IsZombie()) {
164     printf("The %s file does not exist", base.Data());
165     exit(0);
166     }
167    
168    
169     /*
170     * process Mcmd
171     */
172     long int recEntries;
173     tr = (TTree*)rootFile->Get("Mcmd");
174     tr->SetBranchAddress("Mcmd", &mcmdev);
175     nevents = tr->GetEntries();
176    
177     bool timeFound = FALSE;
178 kusanagi 1.4 while (i < nevents) {
179 kusanagi 1.1 tr->GetEntry(i);
180     recEntries = mcmdev->Records->GetEntries();
181 kusanagi 1.4 while (j < recEntries){
182 kusanagi 1.1 mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);
183 kusanagi 1.5 //mcmddata = mcmdrc->McmdData;
184     //printf(" timesync TimeSync %i \n", (unsigned int)mcmddata->At(0));
185     //It is a TimeSync?
186     if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE0)){
187     mcmddata = mcmdrc->McmdData;
188     timesync = (((unsigned int)mcmddata->At(0)<<24)&0xFF000000) + (((unsigned int)mcmddata->At(1)<<16)&0x00FF0000) + (((unsigned int)mcmddata->At(2)<<8)&0x0000FF00) + (((unsigned int)mcmddata->At(3))&0x000000FF);
189     timesync = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
190     timeFound = TRUE;
191     //printf(" timesync TimeSync %i \n", timesync);
192     }
193    
194     //It is an Inclination Mcmd?
195     if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE2)){
196     mcmddata = mcmdrc->McmdData;
197     timesync = (((mcmddata->At(0) << 24) & 0xFF000000) + ((mcmddata->At(1) << 16) & 0x00FF0000) + ((mcmddata->At(2) << 8) & 0x0000FF00) + (mcmddata->At(3) & 0x000000FF))/128.0;
198     timesync = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
199     timeFound = TRUE;
200     //printf(" timesync Inclination %16.8f \n", timesync);
201 kusanagi 1.1 }
202 kusanagi 1.5
203     if (timeFound) break;
204     j++;
205 kusanagi 1.1 }
206 kusanagi 1.5 if (timeFound) break;
207 kusanagi 1.1 i++;
208     }
209    
210 kusanagi 1.4 if (!timeFound) {
211     printf("No timesync info have been found in the file %s", base.Data());
212     exit(0);
213     }
214 kusanagi 1.5
215     //Get the Julian date of the TLE Epoch
216     cJulian offTLE = cJulian(((int)tle1.getField(cTle::FLD_EPOCHYEAR) + 2000), tle1.getField(cTle::FLD_EPOCHDAY));
217    
218     //Get the Julian date of the Resours offset
219     TDatime offRes = TDatime(offDate, offTime);
220     cJulian offResours = cJulian(offRes.GetYear(), offRes.GetMonth(), offRes.GetDay(), offRes.GetHour(), offRes.GetMinute(), offRes.GetSecond());
221     //Add to the Resours Offset the seconds differece beetwen the timesynch from MCMD and the relative OBT
222     offResours.addSec(timesync);
223    
224     //Get the MINUTES past since the TLE Epoch
225     offsetTime = (offResours.getDate() - offTLE.getDate()) * 24.0 * 60.0;
226 kusanagi 1.1 tr = (TTree*)rootFile->Get("Physics");
227     TBranch *headBr = tr->GetBranch("Header");
228     tr->SetBranchAddress("Header", &eh);
229 kusanagi 1.6
230     /********** Anticounter **************/
231     pamela::anticounter::AnticounterEvent *antiev = 0;
232     TH2F *antiRate40 = new TH2F("rate40", base, 360, -180, 180, 161, -80.5, 80.5);
233     TH2F *antiRate30 = new TH2F("rate30", base, 360, -180, 180, 161, -80.5, 80.5);
234     tr->SetBranchAddress("Anticounter", &antiev);
235     Int_t oldAntiRate40 = 0;
236     Int_t newAntiRate40 = 0;
237     Int_t oldAntiRate30 = 0;
238     Int_t newAntiRate30 = 0;
239     Int_t rateAnti40 = 0;
240     Int_t rateAnti30 = 0;
241     /********** Anticounter **************/
242    
243     /********** Trigger **************/
244     pamela::trigger::TriggerEvent *trigger = 0;
245     TH2F *trigAnd = new TH2F("trigAnd", base, 360, -180, 180, 161, -80.5, 80.5); // (S11+S12)*(S21+S22)*(S31+S32)
246     TH2F *trigS1 = new TH2F("trigS1", base, 360, -180, 180, 161, -80.5, 80.5); //(S11+S12)
247     TH2F *trigS111A = new TH2F("trigS111A", base, 360, -180, 180, 161, -80.5, 80.5); //(S111A)
248     tr->SetBranchAddress("Trigger", &trigger);
249     /********** Trigger **************/
250    
251 kusanagi 1.1 nevents = tr->GetEntries();
252     //Fill variables from root-ple
253     for (i = 0; i < nevents; i++){
254 kusanagi 1.5 //for (i = 0; i < 1000; i++){
255 kusanagi 1.1 tr->GetEntry(i);
256     ph = eh->GetPscuHeader();
257 kusanagi 1.5 absTime = offsetTime + (ph->GetOrbitalTime()*(1./60000.));
258 kusanagi 1.1 orbit.getPosition(absTime, &eci);
259     coo = eci.toGeo();
260 kusanagi 1.6 /********** Anticounter **************/
261     oldAntiRate40 = newAntiRate40;
262     oldAntiRate30 = newAntiRate30;
263     newAntiRate40 = antiev->counters[0][6];
264     newAntiRate30 = antiev->counters[0][10];
265     if (newAntiRate40 < oldAntiRate40){
266     newAntiRate40 += 0xFFFF;
267     }
268     if (newAntiRate30 < oldAntiRate30){
269     newAntiRate30 += 0xFFFF;
270     }
271     rateAnti40 = newAntiRate40 - oldAntiRate40;
272     rateAnti30 = newAntiRate30 - oldAntiRate30;
273     /********** Anticounter **************/
274    
275     if (coo.m_Lon > PI) {
276     rate->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat));
277     antiRate30->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), rateAnti30);
278     antiRate40->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), rateAnti40);
279     trigAnd->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), (1/4.)*trigger->trigrate[1]);
280     trigS1->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), 1.*trigger->pmtpl[0]);
281     trigS111A->Fill(rad2deg(coo.m_Lon - 2*PI) , rad2deg(coo.m_Lat), 1.*trigger->pmtcount1[0]);
282     } else {
283     rate->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat));
284     antiRate30->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), rateAnti30);
285     antiRate40->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), rateAnti40);
286     trigAnd->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), (1/4.)*trigger->trigrate[1]);
287     trigS1->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), 1.*trigger->pmtpl[0]);
288     trigS111A->Fill(rad2deg(coo.m_Lon) , rad2deg(coo.m_Lat), 1.*trigger->pmtcount1[0]);
289     //printf(" Time: %16.8f Lat: %16.8f Long: %16.8f Alt: %16.8f\n", absTime, rad2deg(coo.m_Lat), rad2deg(coo.m_Lon), coo.m_Alt);
290     }
291 kusanagi 1.1 }
292     double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000;
293    
294     TMarker *tma=NULL;
295 kusanagi 1.2 TImage *tImage=TImage::Open(mapFile);
296 kusanagi 1.1 int width=(int)(tImage->GetWidth()*0.80);
297     int height=(int)(tImage->GetHeight()*0.80);
298     InitStyle();
299     TCanvas *c1 = new TCanvas("c1","rate/orbit",-width, height); // - : removes the menu bars
300 kusanagi 1.6
301 kusanagi 1.1 TH1F *hframe=NULL;
302     hframe=gPad->DrawFrame(-180,-90,180,90);
303 kusanagi 1.6
304 kusanagi 1.1 oss.str("");
305     oss << filename << " - Event Rate (Hz)";
306     hframe->SetTitle(oss.str().c_str());
307     hframe->SetXTitle("Longitude (deg)");
308     hframe->SetYTitle("Latitude (deg)");
309    
310 kusanagi 1.6 /********** Anticounter **************/
311     TCanvas *anti30 = new TCanvas("anti30","rate/orbit",-width, height); // - : removes the menu bars
312     TH1F *hframe30=NULL;
313     hframe30=gPad->DrawFrame(-180,-90,180,90);
314    
315     oss.str("");
316     oss << filename << " - CAS3 rate/degree";
317     hframe30->SetTitle(oss.str().c_str());
318     hframe30->SetXTitle("Longitude (deg)");
319     hframe30->SetYTitle("Latitude (deg)");
320    
321     TCanvas *anti40 = new TCanvas("anti40","rate/orbit",-width, height); // - : removes the menu bars
322     TH1F *hframe40=NULL;
323     hframe40=gPad->DrawFrame(-180,-90,180,90);
324    
325     oss.str("");
326     oss << filename << " - CAS4 rate/degree";
327     hframe40->SetTitle(oss.str().c_str());
328     hframe40->SetXTitle("Longitude (deg)");
329     hframe40->SetYTitle("Latitude (deg)");
330     /********** Anticounter **************/
331    
332     /********** Trigger **************/
333     TCanvas *trigAndC = new TCanvas("AND","rate/orbit",-width, height); // - : removes the menu bars
334     TH1F *trigAndF=NULL;
335     trigAndF=gPad->DrawFrame(-180,-90,180,90);
336    
337     oss.str("");
338     oss << filename << " - (S11+S12)*(S21+S22)*(S31+S32) rate/degree";
339     trigAndF->SetTitle(oss.str().c_str());
340     trigAndF->SetXTitle("Longitude (deg)");
341     trigAndF->SetYTitle("Latitude (deg)");
342    
343     TCanvas *trigS1C = new TCanvas("trigS1","rate/orbit",-width, height); // - : removes the menu bars
344     TH1F *trigS1F=NULL;
345     trigS1F=gPad->DrawFrame(-180,-90,180,90);
346 kusanagi 1.1
347 kusanagi 1.6 oss.str("");
348     oss << filename << " - (S11+S12) rate/degree";
349     trigS1F->SetTitle(oss.str().c_str());
350     trigS1F->SetXTitle("Longitude (deg)");
351     trigS1F->SetYTitle("Latitude (deg)");
352    
353     TCanvas *trigS111AC = new TCanvas("trigS111A","rate/orbit",-width, height); // - : removes the menu bars
354     TH1F *trigS111AF=NULL;
355     trigS111AF=gPad->DrawFrame(-180,-90,180,90);
356 kusanagi 1.1
357 kusanagi 1.6 oss.str("");
358     oss << filename << " - S111A rate/degree";
359     trigS111AF->SetTitle(oss.str().c_str());
360     trigS111AF->SetXTitle("Longitude (deg)");
361     trigS111AF->SetYTitle("Latitude (deg)");
362     /********** Trigger **************/
363    
364     c1->cd();
365 kusanagi 1.1 TPad *p2 = new TPad("p2", "p2", 0.10, 0.04, 0.983, 1);
366 kusanagi 1.6
367 kusanagi 1.1 p2->cd();
368     TPaletteAxis *tpa=NULL;
369     TH2F *forpal=new TH2F("forpal","",2,0,2,2,0,2);
370     forpal->SetAxisColor(kBlack); //Delete the stat box
371     forpal->SetStats(0); //Delete the stat box
372     forpal->SetMinimum(0.1);
373 kusanagi 1.6 forpal->SetMaximum(rate->GetMaximum());
374     p2->SetLogy();
375 kusanagi 1.1 forpal->SetBinContent(5,1); // just to initialize the histo
376     forpal->SetContour(50);
377     TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette"));
378     tpa=new TPaletteAxis(184,-90,195,90,forpal);
379     tpa->SetLabelColor(kWhite);
380     forpal->Draw("colz");
381 kusanagi 1.6
382     c1->cd();
383     p2->DrawClone();
384    
385     /********** Anticounter **************/
386     Double_t stepPP30 = (antiRate30->GetMaximum() / 50.0);
387     forpal->SetMaximum(antiRate30->GetMaximum());
388     anti30->cd();
389     p2->DrawClone();
390    
391     Double_t stepPP40 = (antiRate40->GetMaximum() / 50.0);
392     forpal->SetMaximum(antiRate40->GetMaximum());
393     anti40->cd();
394     p2->DrawClone();
395     /********** Anticounter **************/
396    
397     /********** Trigger **************/
398     Double_t trigAndPP = (trigAnd->GetMaximum() / 50.0);
399     forpal->SetMaximum(trigAnd->GetMaximum());
400     trigAndC->cd();
401     p2->DrawClone();
402    
403     Double_t trigS1PP = (trigS1->GetMaximum() / 50.0);
404     forpal->SetMaximum(trigS1->GetMaximum());
405     trigS1C->cd();
406     p2->DrawClone();
407    
408     Double_t trigS111APP = (trigS111A->GetMaximum() / 50.0);
409     forpal->SetMaximum(trigS111A->GetMaximum());
410     trigS111AC->cd();
411     p2->DrawClone();
412     /********** Trigger **************/
413 kusanagi 1.1
414     c1->cd();
415     TPad *p1 = new TPad("p1", "p1", 0.10,0.10,0.90,0.92);
416     p1->Draw();
417 kusanagi 1.6
418 kusanagi 1.1 p1->cd();
419     tImage->Draw("X");
420 kusanagi 1.6
421     /********** Anticounter **************/
422     anti30->cd();
423     p1->DrawClone();
424     anti40->cd();
425     p1->DrawClone();
426     /********** Anticounter **************/
427    
428     /********** Trigger **************/
429     trigAndC->cd();
430     p1->DrawClone();
431     trigS1C->cd();
432     p1->DrawClone();
433     trigS111AC->cd();
434     p1->DrawClone();
435     /********** Trigger **************/
436    
437 kusanagi 1.1 c1->cd();
438     gPad->RedrawAxis();
439     c1->Update();
440    
441 kusanagi 1.6 /********** Anticounter **************/
442     anti30->cd();
443     gPad->RedrawAxis();
444     anti30->Update();
445     anti40->cd();
446     gPad->RedrawAxis();
447     anti40->Update();
448     /********** Anticounter **************/
449    
450     /********** Trigger **************/
451     trigAndC->cd();
452     gPad->RedrawAxis();
453     trigAndC->Update();
454    
455     trigS1C->cd();
456     gPad->RedrawAxis();
457     trigS1C->Update();
458    
459     trigS111AC->cd();
460     gPad->RedrawAxis();
461     trigS111AC->Update();
462     /********** Trigger **************/
463    
464 kusanagi 1.1 c1->cd();
465     tma=new TMarker();
466 kusanagi 1.6 tma->SetMarkerStyle(21);
467     tma->SetMarkerSize(0.5);
468     Double_t color = 0.0;
469     Double_t freq = 0.0;
470     Double_t stepPP = (rate->GetMaximum() / 50.0);
471     //Double_t cutoffFreq = stepPP/1000.0;
472     //Double_t alpha = log10(cutoffFreq/stepPP) * (-1.0); //Obviously is 3.0
473    
474     /********** Anticounter **************/
475     Double_t color30 = 0.0;
476     Double_t freq30 = 0.0;
477     //Double_t cutoffFreq30 = stepPP30/1000.0;
478     //Double_t alpha30 = log10(cutoffFreq30/stepPP30) * (-1.0); //Obviously is 3.0
479    
480     Double_t color40 = 0.0;
481     Double_t freq40 = 0.0;
482     //Double_t cutoffFreq40 = stepPP40/1000.0;
483     //Double_t alpha40 = log10(cutoffFreq40/stepPP40) * (-1.0); //Obviously is 3.0
484     /********** Anticounter **************/
485    
486     /********** Trigger **************/
487     Double_t trigAndCol = 0.0;
488     Double_t trigAndFreq = 0.0;
489    
490     Double_t trigS1Col = 0.0;
491     Double_t trigS1Freq = 0.0;
492    
493     Double_t trigS111ACol = 0.0;
494     Double_t trigS111AFreq = 0.0;
495     /********** Trigger **************/
496    
497     for (Int_t kk = 0; kk < 360; kk++){
498     for (Int_t jj = 0; jj < 161; jj++){
499 kusanagi 1.1 freq = rate->GetBinContent(kk, jj);
500 kusanagi 1.6
501     /********** Anticounter **************/
502     freq30 = antiRate30->GetBinContent(kk, jj);
503     freq40 = antiRate40->GetBinContent(kk, jj);
504     /********** Anticounter **************/
505    
506     /********** Trigger **************/
507     trigAndFreq = trigAnd->GetBinContent(kk, jj);
508     trigS1Freq = trigS1->GetBinContent(kk, jj);
509     trigS111AFreq = trigS111A->GetBinContent(kk, jj);
510     /********** Trigger **************/
511    
512     //printf(" Freq: %16.8f \n", freq);
513     posx=kk - 180;
514     posy= jj - 80.5;
515     if (freq > 0) {
516     c1->cd();
517 kusanagi 1.1 // color: palette colors from 51 to 100 ie 50 levels
518 kusanagi 1.6 //color = ((log10(freq/stepPP))/alpha)*49.0 + 100.0; // step defined by palette
519     color = 50.0 + (freq / stepPP); // step defined by palette
520     //printf(" Color: %16.8f \n", color);
521     //if (freq < stepPP) color = 51;
522     //if (color < -100.0) color = 100.0; // just in case if max rate is not max...
523     tma->SetMarkerColor((int)color);
524     tma->DrawMarker(posx,posy);
525     }
526    
527     /********** Anticounter **************/
528     if (freq30 > 0) {
529     anti30->cd();
530     color30 = 50.0 + (freq30 / stepPP30); // step defined by palette
531     //printf(" Color: %16.8f \n", color30);
532     //color30 = ((log10(freq30/stepPP30))/alpha30)*49.0 + 100.0; // step defined by palette
533     //printf(" Color: %16.8f Step: %16.8f Freq: %16.8f \n", color, stepPP, freq);
534     //if (freq30 < cutoffFreq30) color30 = 51;
535     //if (color30 > 100.0) color30 = 100.0; // just in case if max rate is not max...
536     tma->SetMarkerColor((int)color30);
537     tma->DrawMarker(posx,posy);
538     }
539    
540     if (freq40 > 0) {
541     anti40->cd();
542     color40 = 50.0 + (freq40 / stepPP40); // step defined by palette
543     //color40 = ((log10(freq40/stepPP40))/alpha40)*49.0 + 100.0; // step defined by palette
544     //printf(" Color: %16.8f Step: %16.8f Freq: %16.8f \n", color, stepPP, freq);
545     //if (freq40 < cutoffFreq40) color40 = 51;
546     //if (color40 > 100.0) color40 = 100.0; // just in case if max rate is not max...
547     tma->SetMarkerColor((int)color40);
548     tma->DrawMarker(posx,posy);
549     }
550     /********** Anticounter **************/
551    
552     /********** Trigger **************/
553     if (trigAndFreq > 0) {
554     trigAndC->cd();
555     trigAndCol = 50.0 + (trigAndFreq / trigAndPP); // step defined by palette
556     tma->SetMarkerColor((int)trigAndCol);
557 kusanagi 1.1 tma->DrawMarker(posx,posy);
558     }
559 kusanagi 1.6
560     if (trigS1Freq > 0) {
561     trigS1C->cd();
562     trigS1Col = 50.0 + (trigS1Freq / trigS1PP); // step defined by palette
563     tma->SetMarkerColor((int)trigS1Col);
564     tma->DrawMarker(posx,posy);
565     }
566    
567     if (trigS111AFreq > 0) {
568     trigS111AC->cd();
569     trigS111ACol = 50.0 + (trigS111AFreq / trigS111APP); // step defined by palette
570     tma->SetMarkerColor((int)trigS111ACol);
571     tma->DrawMarker(posx,posy);
572     }
573     /********** Anticounter **************/
574     }
575     }
576    
577 kusanagi 1.1 oss.str("");
578     if (outDirectory == "") {
579     oss << filename.Data() << "_OrbitRate." << format.Data();
580     } else {
581     oss << outDirectory.Data() << filename.Data() << "_OrbitRate." << format.Data();
582     }
583 kusanagi 1.6 c1->SaveAs(oss.str().c_str());
584    
585     oss.str("");
586     if (outDirectory == "") {
587     oss << filename.Data() << "_CAS3." << format.Data();
588     } else {
589     oss << outDirectory.Data() << filename.Data() << "_CAS3." << format.Data();
590     }
591     anti30->SaveAs(oss.str().c_str());
592    
593     oss.str("");
594     if (outDirectory == "") {
595     oss << filename.Data() << "_CAS4." << format.Data();
596     } else {
597     oss << outDirectory.Data() << filename.Data() << "_CAS4." << format.Data();
598     }
599     anti40->SaveAs(oss.str().c_str());
600    
601     oss.str("");
602     if (outDirectory == "") {
603     oss << filename.Data() << "_TrigAND." << format.Data();
604     } else {
605     oss << outDirectory.Data() << filename.Data() << "_TrigAND." << format.Data();
606     }
607     trigAndC->SaveAs(oss.str().c_str());
608    
609     oss.str("");
610     if (outDirectory == "") {
611     oss << filename.Data() << "_TrigS1." << format.Data();
612     } else {
613     oss << outDirectory.Data() << filename.Data() << "_TrigS1." << format.Data();
614     }
615     trigS1C->SaveAs(oss.str().c_str());
616    
617     oss.str("");
618     if (outDirectory == "") {
619     oss << filename.Data() << "_trigS111A." << format.Data();
620     } else {
621     oss << outDirectory.Data() << filename.Data() << "_trigS111A." << format.Data();
622     }
623     trigS111AC->SaveAs(oss.str().c_str());
624    
625 kusanagi 1.1 rootFile->Close();
626     }
627    
628 kusanagi 1.6
629    
630    
631 kusanagi 1.1 int main(int argc, char* argv[]){
632 kusanagi 1.2 TString outDir = "./";
633     TString mapFile = "";
634 kusanagi 1.3 TString tleFile = "";
635 kusanagi 1.6 TString format = "png";
636 kusanagi 1.5 int offDate = 20060614;
637     int offTime = 210000;
638 kusanagi 1.1
639     if (argc < 2){
640 kusanagi 1.2 printf("You have to insert at least the file to analyze and the mapFile \n");
641 kusanagi 1.1 printf("Try '--help' for more information. \n");
642     exit(1);
643     }
644    
645     if (!strcmp(argv[1], "--help")){
646 kusanagi 1.2 printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");
647     printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");
648 kusanagi 1.1 printf( "\t --help Print this help and exit \n");
649 kusanagi 1.3 printf( "\t -tleFile[path] Path where to find the tle infos [default dummyOrbit] \n");
650     printf( "\t\tThe tle file have to satisfy a 3-line structure like (this is the included dummyOrbit)\n ");
651     printf( "\t\t\tGP4 Test\n");
652     printf( "\t\t\t1 25544U 98067A 06117.32388940 .00009459 00000-0 64427-4 0 8131\n");
653     printf( "\t\t\t2 25544 51.6388 89.2928 0009043 155.3293 354.6512 15.75673697425143\n");
654 kusanagi 1.1 printf( "\t -outDir[path] Path where to put the output [default ~/tmp] \n");
655 kusanagi 1.6 printf( "\t -format[jpg|gif|ps] Format for output files [default 'png'] \n");
656 kusanagi 1.5 printf( "\t -offDate Date of resetting of the Resource counter [format YYMMDD default 20060614] \n");
657     printf( "\t -offTime Time of resetting of the Resource counter [format HHMMSS default 210000] \n");
658 kusanagi 1.1 exit(1);
659     }
660    
661     for (int i = 2; i < argc; i++){
662     if (!strcmp(argv[i], "-outDir")){
663     if (++i >= argc){
664     printf( "-outDir needs arguments. \n");
665     printf( "Try '--help' for more information. \n");
666     exit(1);
667     } else {
668     outDir = argv[i];
669     continue;
670     }
671     }
672    
673 kusanagi 1.3 if (!strcmp(argv[i], "-tle")){
674     if (++i >= argc){
675     printf( "-tle needs arguments. \n");
676     printf( "Try '--help' for more information. \n");
677     exit(1);
678     } else {
679     tleFile = argv[i];
680     continue;
681     }
682     }
683    
684 kusanagi 1.5 if (!strcmp(argv[i], "-offTime")){
685     if (++i >= argc){
686     printf( "-offTime needs arguments. \n");
687     printf( "Try '--help' for more information. \n");
688     exit(1);
689     }
690     else{
691     offTime = atol(argv[i]);
692     continue;
693     }
694     }
695    
696     if (!strcmp(argv[i], "-offDate")){
697     if (++i >= argc){
698     printf( "-offDate needs arguments. \n");
699     printf( "Try '--help' for more information. \n");
700     exit(1);
701     }
702     else{
703     offDate = atol(argv[i]);
704     continue;
705     }
706     }
707    
708 kusanagi 1.2 if (!strcmp(argv[i], "-map")){
709     if (++i >= argc){
710     printf( "-map needs arguments. \n");
711     printf( "Try '--help' for more information. \n");
712     exit(1);
713     } else {
714     mapFile = argv[i];
715     continue;
716     }
717     }
718    
719     if (!strcmp(argv[i], "-format")) {
720 kusanagi 1.1 if (++i >= argc){
721     printf( "-format needs arguments. \n");
722     printf( "Try '--help' for more information. \n");
723     exit(1);
724     } else {
725     format = argv[i];
726     continue;
727     }
728     }
729     }
730 kusanagi 1.2 if (mapFile != ""){
731 kusanagi 1.5 Rate(argv[1], outDir, format, mapFile, tleFile, offDate, offTime);
732 kusanagi 1.2 } else {
733     printf("You have to insert at least the file to analyze and the mapFile \n");
734     printf("Try '--help' for more information. \n");
735     }
736     }
737    
738 kusanagi 1.1
739    
740    

  ViewVC Help
Powered by ViewVC 1.1.23