quicklook/dataToXML/OrbitalRate.cpp

/**
* Rate
* author  Nagni
* version 1.0 - 27 April 2006
*
* Description: .
*
* Parameters:
*       base - the path where to find the PAMELA unpacked root file.
*       outDirectory - the path where to put the output file.
*       orbPath - the path where to find an orb format for the output.
*
* version 1.0 - 27 April 2006
* First implementation
*/

#include <mcmd/McmdEvent.h>
#include <mcmd/McmdRecord.h>
#include <EventHeader.h>
#include <PscuHeader.h>
#include <TTree.h>
#include "cTle.h"
#include "cEci.h"
#include "cOrbit.h"
#include "TH2F.h"
#include "TFrame.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TASImage.h"
#include "TMarker.h"

#include "TString.h"
#include "TObjString.h"
#include "TStyle.h"
#include "TPaletteAxis.h"
#include "TROOT.h"

#define TRUE  1
#define FALSE 0
/**
 * 
 * @param base 
 * @param outDirectory 
 * @param xslPath 
 */

void InitStyle() {
  gROOT->SetStyle("Plain");

  TStyle *myStyle[2], *tempo;
  myStyle[0]=new TStyle("StyleWhite", "white");
  myStyle[1]=new TStyle("StyleBlack", "black");

  tempo=gStyle;
  Int_t linecol, bkgndcol, histcol;
 
  for(Int_t style=0; style<2; style++) {

    linecol=kWhite*style+kBlack*(1-style);
    bkgndcol=kBlack*style+kWhite*(1-style);
    histcol=kYellow*style+kBlack*(1-style); // was 95

    myStyle[style]->Copy(*tempo);

    myStyle[style]->SetCanvasBorderMode(0);
    myStyle[style]->SetCanvasBorderSize(1);
    myStyle[style]->SetFrameBorderSize(1);
    myStyle[style]->SetFrameBorderMode(0);
    myStyle[style]->SetPadBorderSize(1);
    myStyle[style]->SetStatBorderSize(1);
    myStyle[style]->SetTitleBorderSize(1);
    myStyle[style]->SetPadBorderMode(0);
    myStyle[style]->SetPalette(1,0);
    myStyle[style]->SetPaperSize(20,27);        
    myStyle[style]->SetFuncColor(kRed);
    myStyle[style]->SetFuncWidth(1);
    myStyle[style]->SetLineScalePS(1); 
    myStyle[style]->SetCanvasColor(bkgndcol);
    myStyle[style]->SetAxisColor(linecol,"XYZ");
    myStyle[style]->SetFrameFillColor(bkgndcol);
    myStyle[style]->SetFrameLineColor(linecol);
    myStyle[style]->SetLabelColor(linecol,"XYZ");
    myStyle[style]->SetPadColor(bkgndcol);
    myStyle[style]->SetStatColor(bkgndcol);
    myStyle[style]->SetStatTextColor(linecol);
    myStyle[style]->SetTitleColor(linecol,"XYZ");
    myStyle[style]->SetTitleFillColor(bkgndcol);
    myStyle[style]->SetTitleTextColor(linecol);
    myStyle[style]->SetLineColor(linecol);
    myStyle[style]->SetMarkerColor(histcol);
    myStyle[style]->SetTextColor(linecol);
    
    myStyle[style]->SetGridColor((style)?13:kBlack);
    myStyle[style]->SetHistFillStyle(1001*(1-style));
    myStyle[style]->SetHistLineColor(histcol);
    myStyle[style]->SetHistFillColor((style)?bkgndcol:kYellow);
  }

  myStyle[1]->cd();

  gROOT->ForceStyle();

}

void Rate(TString base, TString outDirectory = "", TString format = "jpg"){
        TTree                  *tr         = 0;
        pamela::McmdEvent      *mcmdev     = 0;
        pamela::McmdRecord     *mcmdrc     = 0;
        TArrayC                *mcmddata   = 0;
        ULong64_t               nevents    = 0;
        ULong64_t               timesync   = 0;
        pamela::EventHeader    *eh         = 0;
        pamela::PscuHeader     *ph         = 0;
        stringstream    oss;
        double offsetTime = 0;
        double absTime;
        UInt_t i = 0;
        UInt_t j = 0;
        
        
        TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
        filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();
        // Test SGP4
        string str1 = "SGP4 Test";
        string str2 = "1 25544U 98067A   06117.32388940  .00009459  00000-0  64427-4 0  8131";
        string str3 = "2 25544  51.6388  89.2928 0009043 155.3293 354.6512 15.75673697425143";

        cTle tle1(str1, str2, str3);
        cOrbit       orbit(tle1);
        cEci         eci;
        cCoordGeo    coo;
        TH2F *rate      = new TH2F("rate", base, 360, -180, 180, 161, -80.5, 80.5);
        TFile *rootFile = new TFile(base);
        if (rootFile->IsZombie()) {
                printf("The %s file does not exist", base.Data());
                exit(0);
        }
        
        
        /* 
         * process Mcmd
         */
        long int recEntries;
        tr = (TTree*)rootFile->Get("Mcmd");
        tr->SetBranchAddress("Mcmd", &mcmdev);
        nevents = tr->GetEntries();
        
        bool timeFound = FALSE;
        while ((i < nevents) | !timeFound){
                tr->GetEntry(i);
                recEntries = mcmdev->Records->GetEntries();
                while ((j < recEntries) | !timeFound){
                        mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);
                        if (mcmdrc->ID1 == 0xE0){
                        mcmddata = mcmdrc->McmdData;
                        timesync = (((ULong64_t)mcmddata->At(0)<<24)&0xFF000000) + 
                                        (((ULong64_t)mcmddata->At(1)<<16)&0x00FF0000) + 
                                        (((ULong64_t)mcmddata->At(2)<<8)&0x0000FF00)  + 
                                        (((ULong64_t)mcmddata->At(3))&0x000000FF);
                        offsetTime = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
                        timeFound = TRUE;
                        }
                j++;
                }
        i++;
        }
        
        tr = (TTree*)rootFile->Get("Physics");
        TBranch *headBr = tr->GetBranch("Header");
        tr->SetBranchAddress("Header", &eh);
        nevents = tr->GetEntries();
        //Fill variables from root-ple
        for (i = 0; i < nevents; i++){
                tr->GetEntry(i);
                ph = eh->GetPscuHeader();
                absTime = ((ph->GetOrbitalTime()*(1./1000.)) + absTime)/60;
                orbit.getPosition(absTime, &eci);
                coo = eci.toGeo();
                rate->Fill(rad2deg(coo.m_Lon), rad2deg(coo.m_Lat));
                /*
                printf("             %16.8f %16.8f %16.8f\n",
                rad2deg(coo.m_Lat),
                rad2deg(coo.m_Lon),
                coo.m_Alt);     
                */
        }
        double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000;

        int ptcnt=0, color=0;
        TMarker *tma=NULL;
        TLine *tli=NULL;
        double step=0;

        TImage *tImage=TImage::Open("Data/terra4.png");
        int width=(int)(tImage->GetWidth()*0.80);
        int height=(int)(tImage->GetHeight()*0.80);
        InitStyle();
        TCanvas *c1 = new TCanvas("c1","rate/orbit",-width, height); // - : removes the menu bars
        
        TH1F *hframe=NULL;
        hframe=gPad->DrawFrame(-180,-90,180,90);
        oss.str("");
        oss << filename << " - Event Rate (Hz)";
        hframe->SetTitle(oss.str().c_str());
        hframe->SetXTitle("Longitude (deg)");
        hframe->SetYTitle("Latitude (deg)");
        c1->cd();
        
        
        TPad *p2 = new TPad("p2", "p2", 0.10, 0.04, 0.983, 1);
        p2->Draw();
        p2->cd();
        TPaletteAxis *tpa=NULL;
        TH2F *forpal=new TH2F("forpal","",2,0,2,2,0,2);
        forpal->SetAxisColor(kBlack); //Delete the stat box
        forpal->SetStats(0); //Delete the stat box
        forpal->SetMinimum(0.1);
        forpal->SetMaximum(200);
        forpal->SetBinContent(5,1);    // just to initialize the histo
        forpal->SetContour(50);
        TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette"));
        step=forpal->GetContourLevel(1)-forpal->GetContourLevel(0);
        tpa=new TPaletteAxis(184,-90,195,90,forpal);
        tpa->SetLabelColor(kWhite);
        forpal->Draw("colz");

        c1->cd();
        TPad *p1 = new TPad("p1", "p1", 0.10,0.10,0.90,0.92);
        p1->Draw();
        p1->cd();
        tImage->Draw("X");
        c1->cd();
        gPad->RedrawAxis();
        c1->Update();
        
  c1->cd();
  ptcnt=0;
  tma=new TMarker();
  tma->SetMarkerStyle(4);
  tli=new TLine();
  tli->SetLineColor(kMagenta);
  Stat_t freq;
      for (Int_t kk = 0; kk < 360; kk++){
        for (Int_t jj = 0; jj < 161; jj++){
                freq = rate->GetBinContent(kk, jj);
                if (freq>0) {
                        posx=(kk - 180); posy= jj - 80.5;
                        // color: palette colors from 51 to 100 ie 50 levels
                        color=51+(int) ((log10((rate==0)?0.1:freq)-log10(0.1))/step); // step defined by palette
                        if (color>100) color=100; // just in case if max rate is not max...
                        tma->SetMarkerColor(color);
                        if (!(posx<0 && oldposx>0) && oldposx!=-1000 && oldposy!=-1000) tli->DrawLine(oldposx,oldposy,posx,posy); 
                        tma->DrawMarker(posx,posy);
                        oldposx=posx;
                        oldposy=posy;   
                }
        }
      }
        oss.str("");
        if (outDirectory == "") {
                oss << filename.Data() << "_OrbitRate." << format.Data();
        } else {
                oss << outDirectory.Data() << filename.Data() << "_OrbitRate." << format.Data();
        }
        c1->SaveAs(oss.str().c_str());
        rootFile->Close();
}

int main(int argc, char* argv[]){
    TString outDir = "./";
    TString format = "jpg";
    
 if (argc < 2){
    printf("You have to insert at least the file to analyze \n");
    printf("Try '--help' for more information. \n");
    exit(1);
  }  

  if (!strcmp(argv[1], "--help")){
        printf( "Usage: LogToXML FILE [OPTION] \n");
        printf( "\t --help                 Print this help and exit \n");
        printf( "\t -outDir[path]          Path where to put the output [default ~/tmp] \n");
        printf( "\t -format[jpg|gif|ps]    Format for output files [default 'jpg'] \n");
        exit(1);
  }

  for (int i = 2; i < argc; i++){
    if (!strcmp(argv[i], "-outDir")){
        if (++i >= argc){
            printf( "-outDir needs arguments. \n");
            printf( "Try '--help' for more information. \n");
            exit(1);
        } else {
            outDir = argv[i];
            continue;
        }
    }

    if (!strcmp(argv[i], "-format")) 
        if (++i >= argc){
            printf( "-format needs arguments. \n");
            printf( "Try '--help' for more information. \n");
            exit(1);
        } else {
                format = argv[i];
                continue;
        }
    }
        Rate(argv[1], outDir, format);
}


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