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

Diff of /quicklook/dataToXML/OrbitalRate.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by kusanagi, Tue May 2 14:32:23 2006 UTC revision 1.5 by kusanagi, Fri Jun 30 13:19:25 2006 UTC
# Line 22  Line 22 
22  #include "cTle.h"  #include "cTle.h"
23  #include "cEci.h"  #include "cEci.h"
24  #include "cOrbit.h"  #include "cOrbit.h"
25    #include "cJulian.h"
26  #include "TH2F.h"  #include "TH2F.h"
27  #include "TFrame.h"  #include "TFrame.h"
28  #include "TGraph.h"  #include "TGraph.h"
29  #include "TCanvas.h"  #include "TCanvas.h"
30  #include "TASImage.h"  #include "TASImage.h"
31  #include "TMarker.h"  #include "TMarker.h"
32    #include <TDatime.h>
33    
34  #include "TString.h"  #include "TString.h"
35  #include "TObjString.h"  #include "TObjString.h"
# Line 35  Line 37 
37  #include "TPaletteAxis.h"  #include "TPaletteAxis.h"
38  #include "TROOT.h"  #include "TROOT.h"
39  #include <sys/stat.h>  #include <sys/stat.h>
40    #include <fstream>
41    
42  #define TRUE  1  #define TRUE  1
43  #define FALSE 0  #define FALSE 0
# Line 104  void InitStyle() { Line 107  void InitStyle() {
107    
108  }  }
109    
110  void Rate(TString base, TString outDirectory = "", TString format = "jpg", TString mapFile = ""){  void Rate(TString base, TString outDirectory = "", TString format = "jpg", TString mapFile = "", TString tleFile = "", int offDate = 20060614, int offTime = 210000){
111          TTree                  *tr         = 0;          TTree                  *tr         = 0;
112          pamela::McmdEvent      *mcmdev     = 0;          pamela::McmdEvent      *mcmdev     = 0;
113          pamela::McmdRecord     *mcmdrc     = 0;          pamela::McmdRecord     *mcmdrc     = 0;
114          TArrayC                *mcmddata   = 0;          pamela::EventHeader    *eh         = 0;
115            pamela::PscuHeader     *ph         = 0;
116            TArrayC                *mcmddata;
117          ULong64_t               nevents    = 0;          ULong64_t               nevents    = 0;
118          ULong64_t               timesync   = 0;          double                  timesync   = 0;
         pamela::EventHeader    *eh         = 0;  
         pamela::PscuHeader     *ph         = 0;  
119          stringstream    oss;          stringstream    oss;
120          double offsetTime = 0;          double offsetTime = 0;
121          double absTime;          double absTime = 0;
122          UInt_t i = 0;          UInt_t i = 0;
123          UInt_t j = 0;          UInt_t j = 0;
124          struct stat buf;          struct stat buf;
# Line 129  void Rate(TString base, TString outDirec Line 132  void Rate(TString base, TString outDirec
132                    
133          TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();          TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
134          filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();          filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();
135          // Test SGP4  
136          string str1 = "SGP4 Test";          string str1 = "RESURS-DK 1";
137          string str2 = "1 25544U 98067A   06117.32388940  .00009459  00000-0  64427-4 0  8131";          string str2 = "1 29228U 06021A   06166.45315890 -.00001271  40192-5  00000+0 0    12";
138          string str3 = "2 25544  51.6388  89.2928 0009043 155.3293 354.6512 15.75673697425143";          string str3 = "2 29228 069.9476 065.6885 0106036 060.6946 300.4810 16.02948390    10";
139            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    
149          cTle tle1(str1, str2, str3);          cTle tle1(str1, str2, str3);
150          cOrbit       orbit(tle1);          cOrbit       orbit(tle1);
# Line 155  void Rate(TString base, TString outDirec Line 167  void Rate(TString base, TString outDirec
167          nevents = tr->GetEntries();          nevents = tr->GetEntries();
168                    
169          bool timeFound = FALSE;          bool timeFound = FALSE;
170          while ((i < nevents) | !timeFound){          while (i < nevents) {
171                  tr->GetEntry(i);                  tr->GetEntry(i);
172                  recEntries = mcmdev->Records->GetEntries();                  recEntries = mcmdev->Records->GetEntries();
173                  while ((j < recEntries) | !timeFound){                  while (j < recEntries){
174                          mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);                          mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);
175                          if (mcmdrc->ID1 == 0xE0){                          //mcmddata = mcmdrc->McmdData;
176                          mcmddata = mcmdrc->McmdData;                          //printf(" timesync TimeSync %i \n", (unsigned int)mcmddata->At(0));
177                          timesync = (((ULong64_t)mcmddata->At(0)<<24)&0xFF000000) +                          //It is a TimeSync?
178                                          (((ULong64_t)mcmddata->At(1)<<16)&0x00FF0000) +                          if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE0)){
179                                          (((ULong64_t)mcmddata->At(2)<<8)&0x0000FF00)  +                                  mcmddata = mcmdrc->McmdData;
180                                          (((ULong64_t)mcmddata->At(3))&0x000000FF);                                  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                          offsetTime = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);                                  timesync = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
182                          timeFound = TRUE;                                  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                          }                          }
194                  j++;                          
195                            if (timeFound) break;
196                            j++;
197                  }                  }
198            if (timeFound) break;
199          i++;          i++;
200          }          }
201                    
202            if (!timeFound) {
203                    printf("No timesync info have been found in the file %s", base.Data());
204                    exit(0);
205            }
206            
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          tr = (TTree*)rootFile->Get("Physics");          tr = (TTree*)rootFile->Get("Physics");
228          TBranch *headBr = tr->GetBranch("Header");          TBranch *headBr = tr->GetBranch("Header");
229          tr->SetBranchAddress("Header", &eh);          tr->SetBranchAddress("Header", &eh);
230          nevents = tr->GetEntries();          nevents = tr->GetEntries();
231          //Fill variables from root-ple          //Fill variables from root-ple
232          for (i = 0; i < nevents; i++){          for (i = 0; i < nevents; i++){
233            //for (i = 0; i < 1000; i++){
234                  tr->GetEntry(i);                  tr->GetEntry(i);
235                  ph = eh->GetPscuHeader();                  ph = eh->GetPscuHeader();
236                  absTime = ((ph->GetOrbitalTime()*(1./1000.)) + absTime)/60;                  absTime = offsetTime + (ph->GetOrbitalTime()*(1./60000.));
237                    //printf(" absTime %16.8f \n", absTime);
238                  orbit.getPosition(absTime, &eci);                  orbit.getPosition(absTime, &eci);
239                  coo = eci.toGeo();                  coo = eci.toGeo();
240                  rate->Fill(rad2deg(coo.m_Lon), rad2deg(coo.m_Lat));                  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);
                 printf("             %16.8f %16.8f %16.8f\n",  
                 rad2deg(coo.m_Lat),  
                 rad2deg(coo.m_Lon),  
                 coo.m_Alt);      
                 */  
242          }          }
243          double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000;          double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000;
244    
# Line 226  void Rate(TString base, TString outDirec Line 273  void Rate(TString base, TString outDirec
273          forpal->SetAxisColor(kBlack); //Delete the stat box          forpal->SetAxisColor(kBlack); //Delete the stat box
274          forpal->SetStats(0); //Delete the stat box          forpal->SetStats(0); //Delete the stat box
275          forpal->SetMinimum(0.1);          forpal->SetMinimum(0.1);
276          forpal->SetMaximum(200);          forpal->SetMaximum(15);
277          forpal->SetBinContent(5,1);    // just to initialize the histo          forpal->SetBinContent(5,1);    // just to initialize the histo
278          forpal->SetContour(50);          forpal->SetContour(50);
279          TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette"));          TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette"));
# Line 280  void Rate(TString base, TString outDirec Line 327  void Rate(TString base, TString outDirec
327  int main(int argc, char* argv[]){  int main(int argc, char* argv[]){
328      TString outDir  = "./";      TString outDir  = "./";
329      TString mapFile = "";      TString mapFile = "";
330        TString tleFile = "";
331      TString format  = "jpg";      TString format  = "jpg";
332        int offDate = 20060614;
333        int offTime = 210000;
334            
335   if (argc < 2){   if (argc < 2){
336      printf("You have to insert at least the file to analyze and the mapFile \n");      printf("You have to insert at least the file to analyze and the mapFile \n");
# Line 292  int main(int argc, char* argv[]){ Line 342  int main(int argc, char* argv[]){
342          printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");          printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");
343          printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");          printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");
344          printf( "\t --help                 Print this help and exit \n");          printf( "\t --help                 Print this help and exit \n");
345            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          printf( "\t -outDir[path]          Path where to put the output [default ~/tmp] \n");          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");          printf( "\t -format[jpg|gif|ps]    Format for output files [default 'jpg'] \n");
352            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          exit(1);          exit(1);
355    }    }
356    
# Line 309  int main(int argc, char* argv[]){ Line 366  int main(int argc, char* argv[]){
366          }          }
367      }      }
368    
369       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        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      if (!strcmp(argv[i], "-map")){      if (!strcmp(argv[i], "-map")){
405          if (++i >= argc){          if (++i >= argc){
406              printf( "-map needs arguments. \n");              printf( "-map needs arguments. \n");
# Line 332  int main(int argc, char* argv[]){ Line 424  int main(int argc, char* argv[]){
424      }      }
425  }  }
426  if (mapFile != ""){  if (mapFile != ""){
427          Rate(argv[1], outDir, format, mapFile);          Rate(argv[1], outDir, format, mapFile, tleFile, offDate, offTime);
428  } else {  } else {
429          printf("You have to insert at least the file to analyze and the mapFile \n");          printf("You have to insert at least the file to analyze and the mapFile \n");
430          printf("Try '--help' for more information. \n");          printf("Try '--help' for more information. \n");
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23