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 "cJulian.h"
#include "TH2F.h"
#include "TFrame.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TASImage.h"
#include "TMarker.h"
#include <TDatime.h>

#include "TString.h"
#include "TObjString.h"
#include "TStyle.h"
#include "TPaletteAxis.h"
#include "TROOT.h"
#include <sys/stat.h>
#include <fstream>

#define TRUE  1
#define FALSE 0
/**
 * 
 * @param base 
 * @param outDirectory 
 * @param xslPath 
 */
using namespace std;

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", TString mapFile = "", TString tleFile = "", int offDate = 20060614, int offTime = 210000){
        TTree                  *tr         = 0;
        pamela::McmdEvent      *mcmdev     = 0;
        pamela::McmdRecord     *mcmdrc     = 0;
        pamela::EventHeader    *eh         = 0;
        pamela::PscuHeader     *ph         = 0;
        TArrayC                *mcmddata;
        ULong64_t               nevents    = 0;
        double                  timesync   = 0;
        stringstream    oss;
        double offsetTime = 0;
        double absTime = 0;
        UInt_t i = 0;
        UInt_t j = 0;
        struct stat buf;
        
        i = stat (mapFile.Data(), &buf );
        // If the file does not exists
        if (i != 0){
                printf("The %s file does not exists.", mapFile.Data());
                exit(0);
        }
        
        TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
        filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();

        string str1 = "RESURS-DK 1";
        string str2 = "1 29228U 06021A   06166.45315890 -.00001271  40192-5  00000+0 0    12";
        string str3 = "2 29228 069.9476 065.6885 0106036 060.6946 300.4810 16.02948390    10";
        if (tleFile != ""){
                fstream fileTle(tleFile.Data(),ios::in);
                if (fileTle.is_open()) {
                        getline (fileTle,str1);
                        getline (fileTle,str2);
                        getline (fileTle,str3);
                }
                fileTle.close();
        }

        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) {
                tr->GetEntry(i);
                recEntries = mcmdev->Records->GetEntries();
                while (j < recEntries){
                        mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);
                        //mcmddata = mcmdrc->McmdData;
                        //printf(" timesync TimeSync %i \n", (unsigned int)mcmddata->At(0));
                        //It is a TimeSync?
                        if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE0)){
                                mcmddata = mcmdrc->McmdData;
                                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);
                                timesync = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
                                timeFound = TRUE;
                                //printf(" timesync TimeSync %i \n", timesync);
                        }
                        
                        //It is an Inclination Mcmd?
                        if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE2)){
                                mcmddata = mcmdrc->McmdData;
                                timesync = (((mcmddata->At(0) << 24) & 0xFF000000) + ((mcmddata->At(1) << 16) & 0x00FF0000) + ((mcmddata->At(2) << 8) & 0x0000FF00) + (mcmddata->At(3) & 0x000000FF))/128.0;
                                timesync = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
                                timeFound = TRUE;
                                //printf(" timesync Inclination %16.8f \n", timesync);
                        }
                        
                        if (timeFound) break;
                        j++;
                }
        if (timeFound) break;
        i++;
        }
        
        if (!timeFound) {
                printf("No timesync info have been found in the file %s", base.Data());
                exit(0);
        }
        
        //printf("%d \n", tle1.getField(cTle::FLD_EPOCHYEAR));
        //printf("%d \n", tle1.getField(cTle::FLD_EPOCHDAY));
        //Get the Julian date of the TLE Epoch
        cJulian offTLE = cJulian(((int)tle1.getField(cTle::FLD_EPOCHYEAR) + 2000), tle1.getField(cTle::FLD_EPOCHDAY));
        //cJulian offTLE = cJulian(2006, 178.79019958);
        
        
        //Get the Julian date of the Resours offset
        TDatime offRes = TDatime(offDate, offTime);
        cJulian offResours = cJulian(offRes.GetYear(), offRes.GetMonth(), offRes.GetDay(), offRes.GetHour(), offRes.GetMinute(), offRes.GetSecond());
        //Add to the Resours Offset the seconds differece beetwen the timesynch from MCMD and the relative OBT
        offResours.addSec(timesync);
        //printf(" diff %16.8f \n", timesync);
        
        //Get the MINUTES past since the TLE Epoch
        offsetTime = (offResours.getDate() - offTLE.getDate()) * 24.0 * 60.0;
        printf(" offSET %16.8f \n", offsetTime);

        //printf(" offResours %16.8f \n", offResours.getDate());
        //printf(" offTLE %16.8f \n", offTLE.getDate());
        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++){
        //for (i = 0; i < 1000; i++){
                tr->GetEntry(i);
                ph = eh->GetPscuHeader();
                absTime = offsetTime + (ph->GetOrbitalTime()*(1./60000.));
                //printf(" absTime %16.8f \n", absTime);
                orbit.getPosition(absTime, &eci);
                coo = eci.toGeo();
                rate->Fill(rad2deg(coo.m_Lon) - 180.0 , rad2deg(coo.m_Lat));
                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);
        }
        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(mapFile);
        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(15);
        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 mapFile = "";
    TString tleFile = "";
    TString format  = "jpg";
    int offDate = 20060614;
    int offTime = 210000;
    
 if (argc < 2){
    printf("You have to insert at least the file to analyze and the mapFile \n");
    printf("Try '--help' for more information. \n");
    exit(1);
  }  

  if (!strcmp(argv[1], "--help")){
        printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");
        printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");
        printf( "\t --help                 Print this help and exit \n");
        printf( "\t -tleFile[path]         Path where to find the tle infos [default dummyOrbit] \n");
        printf( "\t\tThe tle file have to satisfy a 3-line structure like (this is the included dummyOrbit)\n ");
        printf( "\t\t\tGP4 Test\n");
        printf( "\t\t\t1 25544U 98067A   06117.32388940  .00009459  00000-0  64427-4 0  8131\n");
        printf( "\t\t\t2 25544  51.6388  89.2928 0009043 155.3293 354.6512 15.75673697425143\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");
        printf( "\t -offDate               Date of resetting of the Resource counter [format YYMMDD default 20060614] \n");
        printf( "\t -offTime               Time of resetting of the Resource counter [format HHMMSS default 210000] \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], "-tle")){
        if (++i >= argc){
            printf( "-tle needs arguments. \n");
            printf( "Try '--help' for more information. \n");
            exit(1);
        } else {
            tleFile = argv[i];
            continue;
        }
    }

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

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

    if (!strcmp(argv[i], "-map")){
        if (++i >= argc){
            printf( "-map needs arguments. \n");
            printf( "Try '--help' for more information. \n");
            exit(1);
        } else {
            mapFile = 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;
        }
    }
}
if (mapFile != ""){
        Rate(argv[1], outDir, format, mapFile, tleFile, offDate, offTime);
} else {
        printf("You have to insert at least the file to analyze and the mapFile \n");
        printf("Try '--help' for more information. \n");
    }
}


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