/[PAMELA software]/quicklook/SatelliteInclination/src/SatInc.cpp
ViewVC logotype

Diff of /quicklook/SatelliteInclination/src/SatInc.cpp

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

revision 1.1 by cafagna, Thu Feb 8 00:49:36 2007 UTC revision 1.2 by pam-rm2, Tue Apr 17 07:46:48 2007 UTC
# Line 1  Line 1 
1  /**  /**
2   * SaaaatInc   * SaaaatInc
3   * author  Malakhov V.   * author  Malakhov V.
4   * version 1.0 - 05/02/2007   * version 1.0 - 05/02/2007
5   */   */
6    
7  #include "InclinationInfo.h"  #include "InclinationInfo.h"
8    
9  #include <mcmd/McmdEvent.h>  #include <mcmd/McmdEvent.h>
10  #include <mcmd/McmdRecord.h>  #include <mcmd/McmdRecord.h>
11  #include <EventHeader.h>  #include <EventHeader.h>
12  #include <PscuHeader.h>  #include <PscuHeader.h>
13  #include "RunHeaderEvent.h"  #include "RunHeaderEvent.h"
14  #include <TTree.h>  #include <TTree.h>
15  #include "sgp4.h"  #include "sgp4.h"
16  #include "TH2F.h"  #include "TH2F.h"
17    
18  #include "TFrame.h"  #include "TFrame.h"
19  #include "TGraph.h"  #include "TGraph.h"
20  #include "TCanvas.h"  #include "TCanvas.h"
21    
22  #include <TDatime.h>  #include <TDatime.h>
23  #include <TFile.h>  #include <TFile.h>
24    
25  #include <TTimeStamp.h>  #include <TTimeStamp.h>
26  #include "TString.h"  #include "TString.h"
27  #include "TObjString.h"  #include "TObjString.h"
28  #include "TPaletteAxis.h"  #include "TPaletteAxis.h"
29  #include "TROOT.h"  #include "TROOT.h"
30  #include <sys/stat.h>  #include <sys/stat.h>
31  #include <fstream>  #include <fstream>
32  #include <iostream>  #include <iostream>
33    
34  #include "OrbitalRate.h"  #include "OrbitalRate.h"
35    
36  //#include "TMatrixD.h"  //#include "TMatrixD.h"
37    
38    
39  using namespace std;  using namespace std;
40    
41  int main(int argc, char* argv[]){  int main(int argc, char* argv[]){
42    
43    TString *rootFile = NULL;    TString *rootFile = NULL;
44    TString outDir  = "./";    TString outDir  = "./";
45    TString mapFile = "";    TString mapFile = "";
46    TString tleFile = "";    TString tleFile = "";
47    int offDate = 20060928;    int offDate = 20070408;
48      //int offDate = 20060614;      //int offDate = 20060614;
49    int offTime = 210000;    int offTime = 210000;
50    
51    if (argc < 2){    if (argc < 2){
52      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");
53      printf("Try '--help' for more information. \n");      printf("Try '--help' for more information. \n");
54      exit(1);      exit(1);
55    }      }  
56    
57    if (!strcmp(argv[1], "--help")){    if (!strcmp(argv[1], "--help")){
58      //printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");      //printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");
59      //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");
60      printf( "\t --help                 Print this help and exit \n");      printf( "\t --help                 Print this help and exit \n");
61      printf( "\t -tle[File path]         Path where to find the tle infos \n");      printf( "\t -tle[File path]         Path where to find the tle infos \n");
62      printf( "\t\tUse the script retrieve_TLE.sh to create the file.\n ");      printf( "\t\tUse the script retrieve_TLE.sh to create the file.\n ");
63      printf( "\t -outDir[path]          Path where to put the output without last directory.\n");      printf( "\t -outDir[path]          Path where to put the output without last directory.\n");
64      printf( "\t -offDate               Date of resetting of the Resource counter [format YYMMDD (UTC date) default 20060928] \n");      printf( "\t -offDate               Date of resetting of the Resource counter [format YYMMDD (UTC date) default 20060928] \n");
65      printf( "\t -offTime               Time of resetting of the Resource counter [format HHMMSS (UTC date) default 210000] \n");      printf( "\t -offTime               Time of resetting of the Resource counter [format HHMMSS (UTC date) default 210000] \n");
66      exit(1);      exit(1);
67    }    }
68    
69    // Ok, here we should have at least one root file.  We check that    // Ok, here we should have at least one root file.  We check that
70    // the filename contains ".root".    // the filename contains ".root".
71    if(strstr(argv[1], ".root"))    if(strstr(argv[1], ".root"))
72      rootFile = new TString(argv[1]);      rootFile = new TString(argv[1]);
73    else {    else {
74      cerr << "OrbitalRate: no root file." << endl << "See --help" << endl;      cerr << "OrbitalRate: no root file." << endl << "See --help" << endl;
75      exit(EXIT_FAILURE);      exit(EXIT_FAILURE);
76    }    }
77    
78    for (int i = 2; i < argc; i++){    for (int i = 2; i < argc; i++){
79      if (!strcmp(argv[i], "-outDir")){      if (!strcmp(argv[i], "-outDir")){
80        if (++i >= argc){        if (++i >= argc){
81          printf( "-outDir needs arguments. \n");          printf( "-outDir needs arguments. \n");
82          printf( "Try '--help' for more information. \n");          printf( "Try '--help' for more information. \n");
83          exit(1);          exit(1);
84        } else {        } else {
85          outDir = argv[i];          outDir = argv[i];
86          continue;          continue;
87        }        }
88      }      }
89    
90      if (!strcmp(argv[i], "-tle")){      if (!strcmp(argv[i], "-tle")){
91        if (++i >= argc){        if (++i >= argc){
92          printf( "-tle needs arguments. \n");          printf( "-tle needs arguments. \n");
93          printf( "Try '--help' for more information. \n");          printf( "Try '--help' for more information. \n");
94          exit(1);          exit(1);
95        } else {        } else {
96          tleFile = argv[i];          tleFile = argv[i];
97          continue;          continue;
98        }        }
99      }      }
100    
101      if (!strcmp(argv[i], "-offTime")){      if (!strcmp(argv[i], "-offTime")){
102        if (++i >= argc){        if (++i >= argc){
103          printf( "-offTime needs arguments. \n");          printf( "-offTime needs arguments. \n");
104          printf( "Try '--help' for more information. \n");          printf( "Try '--help' for more information. \n");
105          exit(1);          exit(1);
106        }        }
107        else{        else{
108          offTime = atol(argv[i]);          offTime = atol(argv[i]);
109          continue;          continue;
110        }        }
111      }      }
112    
113      if (!strcmp(argv[i], "-offDate")){      if (!strcmp(argv[i], "-offDate")){
114        if (++i >= argc){        if (++i >= argc){
115          printf( "-offDate needs arguments. \n");          printf( "-offDate needs arguments. \n");
116          printf( "Try '--help' for more information. \n");          printf( "Try '--help' for more information. \n");
117          exit(1);          exit(1);
118        }        }
119        else{        else{
120          offDate = atol(argv[i]);          offDate = atol(argv[i]);
121          continue;          continue;
122        }        }
123      }      }
124    
125    }    }
126    
127      Rate(rootFile, outDir, mapFile, tleFile, offDate, offTime);      Rate(rootFile, outDir, mapFile, tleFile, offDate, offTime);
128    
129   }   }
130    
131    
132  void Rate(TString *filename, TString outDirectory = "", TString mapFile = "", TString tleFile = "", int offDate = 20060614, int offTime = 210000)  void Rate(TString *filename, TString outDirectory = "", TString mapFile = "", TString tleFile = "", int offDate = 20060614, int offTime = 210000)
133  {  {
134    // **** Offset to temporarily correct the TDatime bug ****/    // **** Offset to temporarily correct the TDatime bug ****/
135    //  offTime += 10000;    //  offTime += 10000;
136    //********************************************************/    //********************************************************/
137    
138    McmdScan *mcmdReader = new McmdScan();    McmdScan *mcmdReader = new McmdScan();
139    
140    pamela::McmdEvent      *mcmdev     = 0;    pamela::McmdEvent      *mcmdev     = 0;
141    pamela::McmdRecord     *mcmdrc     = 0;    pamela::McmdRecord     *mcmdrc     = 0;
142    pamela::EventHeader    *eh         = 0;    pamela::EventHeader    *eh         = 0;
143    pamela::PscuHeader     *ph       = 0;    pamela::PscuHeader     *ph       = 0;
144    TArrayC                *mcmddata;    TArrayC                *mcmddata;
145        
146    pamela::RunHeaderEvent *reh = new pamela::RunHeaderEvent;    pamela::RunHeaderEvent *reh = new pamela::RunHeaderEvent;
147    pamela::EventHeader    *eH  = new pamela::EventHeader;    pamela::EventHeader    *eH  = new pamela::EventHeader;
148        
149    Quaternions *L_QQ_Q_l = new Quaternions();    Quaternions *L_QQ_Q_l = new Quaternions();
150    InclinationInfo *Plux     = new InclinationInfo();    InclinationInfo *Plux     = new InclinationInfo();
151        
152    float timesync = 0, obt_timesync = 0;    float timesync = 0, obt_timesync = 0;
153        
154    ULong64_t nevents = 0;    ULong64_t nevents = 0;
155    stringstream  oss;    stringstream  oss;
156    
157    Long64_t offsetTime = 0;    Long64_t offsetTime = 0;
158    Long64_t timeElapsedFromTLE = 0;    Long64_t timeElapsedFromTLE = 0;
159    Long64_t deltaTime = 0, oldtimeElapsedFromTLE = 0;    Long64_t deltaTime = 0, oldtimeElapsedFromTLE = 0;
160    bool a_second_is_over;    bool a_second_is_over;
161    
162    // Get a char* to "file" from "/dir1/dir2/.../file.root"    // Get a char* to "file" from "/dir1/dir2/.../file.root"
163    TString basename;    TString basename;
164    basename = ((TObjString*) filename->Tokenize('/')->Last())->GetString(); // we get file.root    basename = ((TObjString*) filename->Tokenize('/')->Last())->GetString(); // we get file.root
165    basename = ((TObjString*)basename.Tokenize('.')->First())->GetString(); // we get file    basename = ((TObjString*)basename.Tokenize('.')->First())->GetString(); // we get file
166        
167    TFile *rootFile = new TFile(*filename);    TFile *rootFile = new TFile(*filename);
168    
169    if (rootFile->IsZombie()) printf("The %s file does not exist",basename.Data());    if (rootFile->IsZombie()) printf("The %s file does not exist",basename.Data());
170    TString fileName = ((TObjString*)basename.Tokenize('/')->Last())->GetString();    TString fileName = ((TObjString*)basename.Tokenize('/')->Last())->GetString();
171    TString filePath = basename;    TString filePath = basename;
172    filePath.ReplaceAll(fileName, "");    filePath.ReplaceAll(fileName, "");
173    filePath.ReplaceAll(".root", "");    filePath.ReplaceAll(".root", "");
174    
175    TTree *tr = (TTree*)rootFile->Get("Mcmd");    TTree *tr = (TTree*)rootFile->Get("Mcmd");
176    nevents = tr->GetEntries();    nevents = tr->GetEntries();
177    tr->SetBranchAddress("Mcmd",&mcmdev);    tr->SetBranchAddress("Mcmd",&mcmdev);
178    tr->SetBranchAddress("Header",&eh);    tr->SetBranchAddress("Header",&eh);
179        
180    ULong_t firstime = 99999999;//ph->GetOrbitalTime();    Long_t firstime = 99999999;//ph->GetOrbitalTime();
181    UInt_t  firstID  = 65535;    UInt_t  firstID  = 65535;
182    ULong_t lastime = 0;//ph->GetOrbitalTime();    Long_t lastime = 0;//ph->GetOrbitalTime();
183    UInt_t lastID = 0;    UInt_t lastID = 0;
184    for(int ii = 0; ii < nevents; ii++){    for(int ii = 0; ii < nevents; ii++){
185      tr->GetEntry(ii);      tr->GetEntry(ii);
186      ph = eh->GetPscuHeader();      ph = eh->GetPscuHeader();
187      Int_t tmpSize = mcmdev->Records->GetEntries();      Int_t tmpSize = mcmdev->Records->GetEntries();
188      for (int qq=0; qq < tmpSize; qq++){      for (int qq=0; qq < tmpSize; qq++){
189          mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(qq);          mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(qq);
190          if (mcmdrc->SeqID <= firstID) firstID = mcmdrc->SeqID;          if (mcmdrc->SeqID <= firstID) firstID = mcmdrc->SeqID;
191          if (mcmdrc->SeqID >= lastID)  lastID  = mcmdrc->SeqID;          if (mcmdrc->SeqID >= lastID)  lastID  = mcmdrc->SeqID;
192          if ((int)mcmdrc->ID1 == 226){          if ((int)mcmdrc->ID1 == 226){
193          L_QQ_Q_l->fill(mcmdrc->McmdData);          L_QQ_Q_l->fill(mcmdrc->McmdData);
194          if (L_QQ_Q_l->time[0] >= lastime) lastime = L_QQ_Q_l->time[0];          if (L_QQ_Q_l->time[0] >= lastime) lastime = L_QQ_Q_l->time[0];
195          if (L_QQ_Q_l->time[0] <= firstime) firstime = L_QQ_Q_l->time[0];          if (L_QQ_Q_l->time[0] <= firstime) firstime = L_QQ_Q_l->time[0];
196          }                        }              
197          };          };
198      }      }
199            
200    TTree *tp = (TTree*)rootFile->Get("RunHeader");    TTree *tp = (TTree*)rootFile->Get("RunHeader");
201    tp->SetBranchAddress("Header", &eH);    tp->SetBranchAddress("Header", &eH);
202    tp->SetBranchAddress("RunHeader", &reh);    tp->SetBranchAddress("RunHeader", &reh);
203    tp->GetEntry(0);    tp->GetEntry(0);
204    ph = eH->GetPscuHeader();    ph = eH->GetPscuHeader();
205    ULong_t TimeSync = reh->LAST_TIME_SYNC_INFO;    ULong_t TimeSync = reh->LAST_TIME_SYNC_INFO;
206          ULong_t ObtSync = reh->OBT_TIME_SYNC;
207    //Get Orientation information    ULong_t DeltaOBT = TimeSync - ObtSync;
208        firstime = firstime - DeltaOBT;
209    TH2F *_L1 = new TH2F("_L1",basename+";OnBoard Time, s; Value of l0",1000,firstime,lastime,1000,-1,1);    lastime = lastime - DeltaOBT;
210    TH2F *_L2 = new TH2F("_L2",basename+";OnBoard Time, s; Value of l1",1000,firstime,lastime,1000,-1,1);      
211    TH2F *_L3 = new TH2F("_L3",basename+";OnBoard Time, s; Value of l2",1000,firstime,lastime,1000,-1,1);    //Get Orientation information
212    TH2F *_L4 = new TH2F("_L4",basename+";OnBoard Time, s; Value of l3",1000,firstime,lastime,1000,-1,1);    
213      TH2F *_L1 = new TH2F("_L1",basename+";OnBoard Time, s; Value of l0",1000,firstime,lastime,1000,-1,1);
214    //TH2F *_A1 = new TH2F("_A1",basename,1000,firstime,lastime,1000,-180,180);    TH2F *_L2 = new TH2F("_L2",basename+";OnBoard Time, s; Value of l1",1000,firstime,lastime,1000,-1,1);
215    //TH2F *_A2 = new TH2F("_A2",basename,1000,firstime,lastime,1000,-180,180);    TH2F *_L3 = new TH2F("_L3",basename+";OnBoard Time, s; Value of l2",1000,firstime,lastime,1000,-1,1);
216    //TH2F *_A3 = new TH2F("_A3",basename,1000,firstime,lastime,1000,-180,180);    TH2F *_L4 = new TH2F("_L4",basename+";OnBoard Time, s; Value of l3",1000,firstime,lastime,1000,-1,1);
217      
218    //TH2F *_secID = new TH2F("_secID",basename,1000,firstime,lastime,1000,firstID,lastID);    //TH2F *_A1 = new TH2F("_A1",basename,1000,firstime,lastime,1000,-180,180);
219        //TH2F *_A2 = new TH2F("_A2",basename,1000,firstime,lastime,1000,-180,180);
220    TH2F *_Tangazh = new TH2F("Pitch",basename+";OnBoard Time, s;Angel of pith, deg",1000,firstime,lastime,1000,-20,20);    //TH2F *_A3 = new TH2F("_A3",basename,1000,firstime,lastime,1000,-180,180);
221    TH2F *_Ryskanie = new TH2F("Yaw",basename+";OnBoard Time, s;Angel of yaw, deg",1000,firstime,lastime,1000,-100,-80);    
222    TH2F *_Kren = new TH2F("Roll",basename+";OnBoard Time, s;Angel of roll, deg",1000,firstime,lastime,1000,-50,50);    //TH2F *_secID = new TH2F("_secID",basename,1000,firstime,lastime,1000,firstID,lastID);
223        
224    // Open the root file.    TH2F *_Tangazh = new TH2F("Pitch",basename+";OnBoard Time, s;Angle of pith, deg",1000,firstime,lastime,1000,-20,20);
225    rootFile = new TFile(filename->Data());    TH2F *_Ryskanie = new TH2F("Yaw",basename+";OnBoard Time, s;Angle of yaw, deg",1000,firstime,lastime,1000,-20,20);
226    if (rootFile->IsZombie()) {    TH2F *_Kren = new TH2F("Roll",basename+";OnBoard Time, s;Angle of roll, deg",1000,firstime,lastime,1000,-50,50);
227      printf("The file %s does not exist\n", (filename->Data()));    
228      exit(EXIT_FAILURE);    // Open the root file.
229    }    rootFile = new TFile(filename->Data());
230      if (rootFile->IsZombie()) {
231    // Look for a timesync in the TFile rootFile.  We also get the obt      printf("The file %s does not exist\n", (filename->Data()));
232    // of the timesync mcmd.      exit(EXIT_FAILURE);
233    bool err;    }
234    err = lookforTimesync(rootFile, &timesync, &obt_timesync);  
235    if(!err) {    // Look for a timesync in the TFile rootFile.  We also get the obt
236      cerr << "Warning!!! No timesync info has been found in the file "    // of the timesync mcmd.
237           << filename->Data() << endl;    bool err;
238      exit(EXIT_FAILURE);    err = lookforTimesync(rootFile, &timesync, &obt_timesync);
239    }    if(!err) {
240        cerr << "Warning!!! No timesync info has been found in the file "
241    //Get the Julian date of the Resours offset           << filename->Data() << endl;
242    TDatime offRes = TDatime(offDate, offTime);      exit(EXIT_FAILURE);
243    // Add to the Resours Offset the timesync.  This is now the date at    }
244    // the moment of the timesync.  
245    ULong_t TH=offRes.Convert();    //Get the Julian date of the Resours offset
246    //cout<<"TH= "<<TH/86400<<"\n";    TDatime offRes = TDatime(offDate, offTime);
247    offRes.Set(offRes.Convert() + (UInt_t) timesync);    // Add to the Resours Offset the timesync.  This is now the date at
248        // the moment of the timesync.
249    // Now I need a pointer to a cTle object.  The class misses a    ULong_t TH=offRes.Convert();
250    // constructor without arguments, so we have to give it a dummy TLE.    //cout<<"TH= "<<TH/86400<<"\n";
251    string str1 = "RESURS-DK 1";    offRes.Set(offRes.Convert() + (UInt_t) timesync);
252    string str2 = "1 29228U 06021A   06170.19643714  .00009962  00000-0  21000-3 0   196";    
253    string str3 = "2 29228 069.9363 054.7893 0167576 127.4359 017.0674 15.31839265   604";    // Now I need a pointer to a cTle object.  The class misses a
254    cTle *tle1 = new cTle(str1, str2, str3);    // constructor without arguments, so we have to give it a dummy TLE.
255      string str1 = "RESURS-DK 1";
256    // If we have to use a TLE file, call getTle().    string str2 = "1 29228U 06021A   06170.19643714  .00009962  00000-0  21000-3 0   196";
257    if (tleFile != "")    string str3 = "2 29228 069.9363 054.7893 0167576 127.4359 017.0674 15.31839265   604";
258      tle1 = getTle(tleFile, offRes);    cTle *tle1 = new cTle(str1, str2, str3);
259    
260    cOrbit       orbit(*tle1);    // If we have to use a TLE file, call getTle().
261    cEci         eci;    if (tleFile != "")
262    cCoordGeo    coo;      tle1 = getTle(tleFile, offRes);
263    
264    // Get the Julian date of the TLE epoch    cOrbit       orbit(*tle1);
265    string datetime = getTleDatetime(tle1);    cEci         eci;
266    TDatime tledate = TDatime(datetime.c_str());    cCoordGeo    coo;
267    
268          // Get the Julian date of the TLE epoch
269    tr = (TTree*)rootFile->Get("Mcmd");    string datetime = getTleDatetime(tle1);
270    nevents = tr->GetEntries();    TDatime tledate = TDatime(datetime.c_str());
271    tr->SetBranchAddress("Mcmd", &mcmdev);  
272    tr->SetBranchAddress("Header", &eh);      
273        tr = (TTree*)rootFile->Get("Mcmd");
274   // oss.str("");    nevents = tr->GetEntries();
275   // oss << basename+".txt";    tr->SetBranchAddress("Mcmd", &mcmdev);
276        tr->SetBranchAddress("Header", &eh);
277   // ofstream outputFile;    
278   // outputFile.open(oss.str().c_str());   // oss.str("");
279       // oss << basename+".txt";
280    Int_t tmpSize;    
281       // ofstream outputFile;
282    for(UInt_t i = 0; i < nevents; i++) //Fill variables from root-ple   // outputFile.open(oss.str().c_str());
283      {    
284        tr->GetEntry(i);    Int_t tmpSize;
285        tmpSize = mcmdev->Records->GetEntries();    
286        ph = eh->GetPscuHeader();        for(UInt_t i = 0; i < nevents; i++) //Fill variables from root-ple
287              {
288        for (Int_t j3 = 0; j3 < tmpSize; j3++){        tr->GetEntry(i);
289          mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);        tmpSize = mcmdev->Records->GetEntries();
290          if ((int)mcmdrc->ID1 == 226){        ph = eh->GetPscuHeader();    
291              L_QQ_Q_l->fill(mcmdrc->McmdData);        
292              Plux->QuaternionstoAngle(*L_QQ_Q_l);        for (Int_t j3 = 0; j3 < tmpSize; j3++){
293              for (Int_t oi = 0; oi < 6; oi++){          mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j3);
294                  //if ((Int_t)L_QQ_Q_l->CodeErrQua[oi]==0){          if ((int)mcmdrc->ID1 == 226){
295                                L_QQ_Q_l->fill(mcmdrc->McmdData);
296                      _L1->Fill(L_QQ_Q_l->time[oi],L_QQ_Q_l->quat[oi][0]);              //Plux->QuaternionstoAngle(*L_QQ_Q_l);
297                      _L2->Fill(L_QQ_Q_l->time[oi],L_QQ_Q_l->quat[oi][1]);              for (Int_t oi = 0; oi < 6; oi++){
298                      _L3->Fill(L_QQ_Q_l->time[oi],L_QQ_Q_l->quat[oi][2]);                  //if ((Int_t)L_QQ_Q_l->CodeErrQua[oi]==0){
299                      _L4->Fill(L_QQ_Q_l->time[oi],L_QQ_Q_l->quat[oi][3]);                      Int_t timeO = L_QQ_Q_l->time[oi] - DeltaOBT;
300                                            _L1->Fill(timeO,L_QQ_Q_l->quat[oi][0]);
301                      //outputFile << " " << L_QQ_Q_l->time[oi] << " " << -L_QQ_Q_l->quat[oi][0] << " " << -L_QQ_Q_l->quat[oi][1] << " " << -L_QQ_Q_l->quat[oi][2] << " " << -L_QQ_Q_l->quat[oi][3];                      _L2->Fill(timeO,L_QQ_Q_l->quat[oi][1]);
302                                            _L3->Fill(timeO,L_QQ_Q_l->quat[oi][2]);
303                      //_A1->Fill(L_QQ_Q_l->time[oi],Plux->A11);                      _L4->Fill(timeO,L_QQ_Q_l->quat[oi][3]);
304                      //_A2->Fill(L_QQ_Q_l->time[oi],Plux->A12);                      
305                      //_A3->Fill(L_QQ_Q_l->time[oi],Plux->A13);                      //outputFile << " " << L_QQ_Q_l->time[oi] << " " << -L_QQ_Q_l->quat[oi][0] << " " << -L_QQ_Q_l->quat[oi][1] << " " << -L_QQ_Q_l->quat[oi][2] << " " << -L_QQ_Q_l->quat[oi][3];
306                                            
307                      //outputFile << " " << Plux->A11 << " " << Plux->A12 << " " << Plux->A13;                      //_A1->Fill(L_QQ_Q_l->time[oi],Plux->A11);
308                                            //_A2->Fill(L_QQ_Q_l->time[oi],Plux->A12);
309                      //_secID->Fill(L_QQ_Q_l->time[oi],mcmdrc->SeqID);                      //_A3->Fill(L_QQ_Q_l->time[oi],Plux->A13);
310                                            
311                      timeElapsedFromTLE = TH - (Long64_t) tledate.Convert()+ L_QQ_Q_l->time[oi];                      //outputFile << " " << Plux->A11 << " " << Plux->A12 << " " << Plux->A13;
312                      orbit.getPosition(((double) timeElapsedFromTLE)/60., &eci);                      
313                      //outputFile<<" "<< eci.getPos().m_z;//eci.getPos().m_x<<" "<<eci.getPos().m_y<<" "<<eci.getPos().m_z<<" "<<eci.getVel().m_x<<" "<<eci.getVel().m_y<<" "<<eci.getVel().m_z;                      //_secID->Fill(L_QQ_Q_l->time[oi],mcmdrc->SeqID);
314                      Plux->TransAngle(eci.getPos().m_x,eci.getPos().m_y,eci.getPos().m_z,eci.getVel().m_x,eci.getVel().m_y,eci.getVel().m_z,Plux->A11,Plux->A12,Plux->A13,L_QQ_Q_l->quat[oi][0],L_QQ_Q_l->quat[oi][1],L_QQ_Q_l->quat[oi][2],L_QQ_Q_l->quat[oi][3]);                      
315                                            timeElapsedFromTLE = TH - (Long64_t) tledate.Convert()+ L_QQ_Q_l->time[oi];
316                      _Tangazh->Fill(L_QQ_Q_l->time[oi],Plux->Tangazh);                      orbit.getPosition(((double) timeElapsedFromTLE)/60., &eci);
317                      _Ryskanie->Fill(L_QQ_Q_l->time[oi],Plux->Ryskanie);                      //outputFile<<" "<< eci.getPos().m_z;//eci.getPos().m_x<<" "<<eci.getPos().m_y<<" "<<eci.getPos().m_z<<" "<<eci.getVel().m_x<<" "<<eci.getVel().m_y<<" "<<eci.getVel().m_z;
318                      _Kren->Fill(L_QQ_Q_l->time[oi],Plux->Kren);                      Plux->TransAngle(eci.getPos().m_x,eci.getPos().m_y,eci.getPos().m_z,eci.getVel().m_x,eci.getVel().m_y,eci.getVel().m_z,L_QQ_Q_l->quat[oi][0],L_QQ_Q_l->quat[oi][1],L_QQ_Q_l->quat[oi][2],L_QQ_Q_l->quat[oi][3]);
319                                            
320                      //outputFile<<" "<<Plux->Tangazh<<"  "<<Plux->Kren<<" "<<Plux->Ryskanie<<"\n";                      _Tangazh->Fill(timeO,Plux->Tangazh);
321                  //}                                  _Ryskanie->Fill(timeO,Plux->Ryskanie);
322              }                      _Kren->Fill(timeO,Plux->Kren);
323          }                      
324        }                      //outputFile<<" "<<Plux->Tangazh<<"  "<<Plux->Kren<<" "<<Plux->Ryskanie<<"\n";
325      }                  //}            
326                }
327      //outputFile.close();          }
328              }
329      TCanvas *c1 = new TCanvas("c1","Quaternions",1200,1600);      }
330      //TCanvas *c2 = new TCanvas("c2","Anglees",1200,1600);  
331      TCanvas *c4 = new TCanvas("c4","AngleesTRK",1200,1600);        //outputFile.close();
332          
333      c1->Divide(1,4);      TCanvas *c1 = new TCanvas("c1","Quaternions",1200,1600);
334      //c2->Divide(1,3);      //TCanvas *c2 = new TCanvas("c2","Anglees",1200,1600);
335      c4->Divide(1,3);      TCanvas *c4 = new TCanvas("c4","AngleesTRK",1200,1600);  
336      
337      c1->cd(1);      c1->Divide(1,4);
338      _L1->SetMarkerColor(kBlack);      //c2->Divide(1,3);
339      _L1->Draw("");      c4->Divide(1,3);
340      
341      c1->cd(2);      c1->cd(1);
342      _L2->SetMarkerColor(kBlue);      _L1->SetMarkerColor(kBlack);
343      _L2->Draw("");      _L1->Draw("");
344          
345      c1->cd(3);      c1->cd(2);
346      _L3->SetMarkerColor(kRed);      _L2->SetMarkerColor(kBlue);
347      _L3->Draw("");      _L2->Draw("");
348          
349      c1->cd(4);      c1->cd(3);
350      _L4->SetMarkerColor(kGreen);      _L3->SetMarkerColor(kRed);
351      _L4->Draw("");      _L3->Draw("");
352        
353  //    c1->cd(5);      c1->cd(4);
354  //    _secID->SetMarkerColor(kYellow);      _L4->SetMarkerColor(kGreen);
355  //    _secID->Draw("");      _L4->Draw("");
356  /*      
357      c2->cd(1);  //    c1->cd(5);
358      _A3->SetMarkerColor(kRed);  //    _secID->SetMarkerColor(kYellow);
359      _A3->Draw("");  //    _secID->Draw("");
360      /*  
361      c2->cd(2);      c2->cd(1);
362      _A2->SetMarkerColor(kGreen);      _A3->SetMarkerColor(kRed);
363      _A2->Draw("");      _A3->Draw("");
364        
365      c2->cd(3);      c2->cd(2);
366      _A1->SetMarkerColor(13);      _A2->SetMarkerColor(kGreen);
367      _A1->Draw("");      _A2->Draw("");
368  */          
369      c4->cd(1);      c2->cd(3);
370      _Tangazh->SetMarkerColor(142);      _A1->SetMarkerColor(13);
371      _Tangazh->Draw("");      _A1->Draw("");
372            */      
373      c4->cd(2);      c4->cd(1);
374      _Ryskanie->SetMarkerColor(226);      _Tangazh->SetMarkerColor(142);
375      _Ryskanie->Draw("");      _Tangazh->Draw("");
376                    
377      c4->cd(3);      c4->cd(2);
378      _Kren->SetMarkerColor(142);      _Ryskanie->SetMarkerColor(226);
379      _Kren->Draw("");      _Ryskanie->Draw("");
380                    
381      oss.str("");      c4->cd(3);
382      oss << outDirectory+basename+"/"+basename << "_Inclinations_Quaternions.png";      _Kren->SetMarkerColor(142);
383      c1->SaveAs(oss.str().c_str());      _Kren->Draw("");
384      oss.str();          
385        oss.str("");
386  /*    oss.str("");      oss << outDirectory+basename+"/"+basename << "_Inclinations_Quaternions.png";
387      oss << basename<< "_Angles.png";      c1->SaveAs(oss.str().c_str());
388      c2->SaveAs(oss.str().c_str());      oss.str();
389      oss.str();  
390    */  /*    oss.str("");
391      oss.str("");      oss << basename<< "_Angles.png";
392      oss <<  outDirectory+basename+"/"+basename << "_Inclinations_Angles.png";      c2->SaveAs(oss.str().c_str());
393      c4->SaveAs(oss.str().c_str());      oss.str();
394      oss.str();    */
395              oss.str("");
396        oss <<  outDirectory+basename+"/"+basename << "_Inclinations_Angles.png";
397        c4->SaveAs(oss.str().c_str());
398    delete _L1;      oss.str();
399    delete _L2;        
400    delete _L3;  
401    delete _L4;  
402        delete _L1;
403   // delete _A1;    delete _L2;
404   // delete _A2;    delete _L3;
405   // delete _A3;    delete _L4;
406        
407    //delete _secID;   // delete _A1;
408       // delete _A2;
409    delete _Tangazh;   // delete _A3;
410    delete _Kren;    
411    delete _Ryskanie;    //delete _secID;
412      
413    rootFile->Close();    delete _Tangazh;
414  }    delete _Kren;
415      delete _Ryskanie;
416  // Get the TLE from tleFile.  The right TLE is that one with the  
417  // closest previous date to offRes, that is the date at the time of    rootFile->Close();
418  // the first timesync found in the root file.  }
419  //  
420  // Warning: you must use a tle file obtained by space-track.org  // Get the TLE from tleFile.  The right TLE is that one with the
421  // querying the database with the RESURS DK-1 id number 29228,  // closest previous date to offRes, that is the date at the time of
422  // selecting the widest timespan, including the satellite name in the  // the first timesync found in the root file.
423  // results.  //
424  cTle *getTle(TString tleFile, TDatime offRes)  // Warning: you must use a tle file obtained by space-track.org
425  {  // querying the database with the RESURS DK-1 id number 29228,
426    Float_t tledatefromfile, tledatefromroot;  // selecting the widest timespan, including the satellite name in the
427    fstream tlefile(tleFile.Data(), ios::in);  // results.
428    vector<cTle*> ctles;  cTle *getTle(TString tleFile, TDatime offRes)
429    vector<cTle*>::iterator iter;  {
430      Float_t tledatefromfile, tledatefromroot;
431      fstream tlefile(tleFile.Data(), ios::in);
432    // Build a vector of *cTle      vector<cTle*> ctles;
433    while(1) {    vector<cTle*>::iterator iter;
434      cTle *tlef;  
435      string str1, str2, str3;  
436      // Build a vector of *cTle  
437      getline(tlefile, str1);    while(1) {
438      if(tlefile.eof()) break;      cTle *tlef;
439        string str1, str2, str3;
440      getline(tlefile, str2);  
441      if(tlefile.eof()) break;      getline(tlefile, str1);
442        if(tlefile.eof()) break;
443      getline(tlefile, str3);  
444      if(tlefile.eof()) break;      getline(tlefile, str2);
445        if(tlefile.eof()) break;
446      // We now have three good lines for a cTle.  
447      tlef = new cTle(str1, str2, str3);      getline(tlefile, str3);
448      ctles.push_back(tlef);      if(tlefile.eof()) break;
449    }  
450        // We now have three good lines for a cTle.
451    // Sort by date      tlef = new cTle(str1, str2, str3);
452    sort(ctles.begin(), ctles.end(), compTLE);      ctles.push_back(tlef);
453      }
454    tledatefromroot = (offRes.GetYear()-2000)*1e3 + (offRes.Convert() - (TDatime(offRes.GetYear(), 1, 1, 0, 0, 0)).Convert())/ (24.*3600.);  
455      // Sort by date
456    for(iter = ctles.begin(); iter != ctles.end(); iter++) {    sort(ctles.begin(), ctles.end(), compTLE);
457      cTle *tle = *iter;  
458      tledatefromroot = (offRes.GetYear()-2000)*1e3 + (offRes.Convert() - (TDatime(offRes.GetYear(), 1, 1, 0, 0, 0)).Convert())/ (24.*3600.);
459      tledatefromfile = getTleJulian(tle);  
460      for(iter = ctles.begin(); iter != ctles.end(); iter++) {
461      if(tledatefromroot > tledatefromfile) {      cTle *tle = *iter;
462        tlefile.close();  
463        cTle *thisTle = tle;      tledatefromfile = getTleJulian(tle);
464        ctles.clear();  
465        if(tledatefromroot > tledatefromfile) {
466        return thisTle;        tlefile.close();
467      }        cTle *thisTle = tle;
468    }        ctles.clear();
469    
470    // File ended withoud founding a TLE with a date after offRes.  We'll use the last aviable date.        return thisTle;
471    cerr << "Warning: using last available TLE in " << tleFile.Data() << ".  Consider updating your tle file." << endl;      }
472      }
473    tlefile.close();  
474    cTle *thisTle = ctles[ctles.size()-1];    // File ended withoud founding a TLE with a date after offRes.  We'll use the last aviable date.
475    ctles.clear();    cerr << "Warning: using last available TLE in " << tleFile.Data() << ".  Consider updating your tle file." << endl;
476    
477    return thisTle;    tlefile.close();
478  }    cTle *thisTle = ctles[ctles.size()-1];
479      ctles.clear();
480    
481  // Return whether the first TLE is older than the second    return thisTle;
482  bool compTLE (cTle *tle1, cTle *tle2)  }
483  {  
484    return getTleJulian(tle1) > getTleJulian(tle2);  
485  }  // Return whether the first TLE is older than the second
486    bool compTLE (cTle *tle1, cTle *tle2)
487    {
488  // Return the date of the tle using the format (year-2000)*1e3 +    return getTleJulian(tle1) > getTleJulian(tle2);
489  // julian day.  e.g. 6364.5 is the 31th Dec 2006 12:00:00.  }
490  // It does *not* return a cJulian date.  
491  float getTleJulian(cTle *tle) {  
492    return tle->getField(cTle::FLD_EPOCHYEAR)*1e3 + tle->getField(cTle::FLD_EPOCHDAY);  // Return the date of the tle using the format (year-2000)*1e3 +
493  }  // julian day.  e.g. 6364.5 is the 31th Dec 2006 12:00:00.
494    // It does *not* return a cJulian date.
495    float getTleJulian(cTle *tle) {
496  // Look for a timesync in the TFile rootFile.  Set timesync and    return tle->getField(cTle::FLD_EPOCHYEAR)*1e3 + tle->getField(cTle::FLD_EPOCHDAY);
497  // obt_timesync.  Returns 1 if timesync is found, 0 otherwise.  }
498  int lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) {  
499    *timesync = -1;  // will be != -1 if found  
500    // Look for a timesync in the TFile rootFile.  Set timesync and
501    ULong64_t             nevents    = 0;  // obt_timesync.  Returns 1 if timesync is found, 0 otherwise.
502    pamela::McmdRecord     *mcmdrc     = 0;  UInt_t lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) {
503    pamela::McmdEvent      *mcmdev     = 0;    *timesync = -1;  // will be != -1 if found
504    TArrayC                *mcmddata;  
505    TTree *tr = (TTree*) rootFile->Get("Mcmd");    ULong64_t             nevents    = 0;
506      pamela::McmdRecord     *mcmdrc     = 0;
507    tr->SetBranchAddress("Mcmd", &mcmdev);    pamela::McmdEvent      *mcmdev     = 0;
508      TArrayC                *mcmddata;
509    nevents = tr->GetEntries();    TTree *tr = (TTree*) rootFile->Get("Mcmd");
510    
511    // Looking for a timesync.  We stop at the first one found.    tr->SetBranchAddress("Mcmd", &mcmdev);
512    long int recEntries;  
513      nevents = tr->GetEntries();
514    for(UInt_t i = 0; i < nevents; i++) {  
515      tr->GetEntry(i);    // Looking for a timesync.  We stop at the first one found.
516      recEntries = mcmdev->Records->GetEntries();    long int recEntries;
517    
518      for(UInt_t j = 0; j < recEntries; j++) {    for(UInt_t i = 0; i < nevents; i++) {
519        mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);      tr->GetEntry(i);
520        recEntries = mcmdev->Records->GetEntries();
521        if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE0)) //Is it a TimeSync?  
522          {      for(UInt_t j = 0; j < recEntries; j++) {
523            mcmddata = mcmdrc->McmdData;        mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);
524            *timesync = (((unsigned int)mcmddata->At(0)<<24)&0xFF000000)  
525              + (((unsigned int)mcmddata->At(1)<<16)&0x00FF0000)        if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE0)) //Is it a TimeSync?
526              + (((unsigned int)mcmddata->At(2)<<8)&0x0000FF00)          {
527              + (((unsigned int)mcmddata->At(3))&0x000000FF);            mcmddata = mcmdrc->McmdData;
528            *obt_timesync = (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);            *timesync = (((unsigned int)mcmddata->At(0)<<24)&0xFF000000)
529                + (((unsigned int)mcmddata->At(1)<<16)&0x00FF0000)
530            goto out;  // a timesync is found              + (((unsigned int)mcmddata->At(2)<<8)&0x0000FF00)
531          }              + (((unsigned int)mcmddata->At(3))&0x000000FF);
532      }            *obt_timesync = (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
533    }  
534              goto out;  // a timesync is found
535   out:          }
536                }
537    if (*timesync == -1)    }
538      return 0;  
539    else   out:
540      return 1;          
541  }    if (*timesync == -1)
542        return 0;
543      else
544  // Return a string like YYYY-MM-DD hh:mm:ss, a datetime format.      return 1;
545  string getTleDatetime(cTle *tle)  }
546  {  
547    int year, mon, day, hh, mm, ss;  
548    double dom; // day of month (is double!)  // Return a string like YYYY-MM-DD hh:mm:ss, a datetime format.
549    stringstream date; // date in datetime format  string getTleDatetime(cTle *tle)
550    {
551    // create a cJulian from the date in tle    int year, mon, day, hh, mm, ss;
552    cJulian jdate = cJulian( 2000 + (int) tle->getField(cTle::FLD_EPOCHYEAR), tle->getField(cTle::FLD_EPOCHDAY));    double dom; // day of month (is double!)
553      stringstream date; // date in datetime format
554    // get year, month, day of month  
555    jdate.getComponent(&year, &mon, &dom);    // create a cJulian from the date in tle
556      cJulian jdate = cJulian( 2000 + (int) tle->getField(cTle::FLD_EPOCHYEAR), tle->getField(cTle::FLD_EPOCHDAY));
557    // build a datetime YYYY-MM-DD hh:mm:ss  
558    date.str("");    // get year, month, day of month
559    day = (int) floor(dom);    jdate.getComponent(&year, &mon, &dom);
560    hh = (int) floor( (dom - day) * 24);  
561    mm = (int) floor( ((dom - day) * 24 - hh) * 60);    // build a datetime YYYY-MM-DD hh:mm:ss
562    ss = (int) floor( ((((dom - day) * 24 - hh) * 60 - mm) * 60));    date.str("");
563    //  ms = (int) floor( (((((dom - day) * 24 - hh) * 60 - mm) * 60) - ss) * 1000);    day = (int) floor(dom);
564      hh = (int) floor( (dom - day) * 24);
565    date << year << "-" << mon << "-" << day << " " << hh << ":" << mm << ":" << ss;    mm = (int) floor( ((dom - day) * 24 - hh) * 60);
566      ss = (int) floor( ((((dom - day) * 24 - hh) * 60 - mm) * 60));
567    return date.str();    //  ms = (int) floor( (((((dom - day) * 24 - hh) * 60 - mm) * 60) - ss) * 1000);
568  }  
569      date << year << "-" << mon << "-" << day << " " << hh << ":" << mm << ":" << ss;
570  //  
571  // Solve the overflow for anticoincidence because this counter is    return date.str();
572  // stored in 2 bytes so counts from 0 to 65535.  }
573  //  
574  // counter is the actual value.  //
575  // oldValue is meant to be the previous value of counter.  // Solve the overflow for anticoincidence because this counter is
576  //  // stored in 2 bytes so counts from 0 to 65535.
577  // Example:  //
578  // for(...) {  // counter is the actual value.
579  //   ...  // oldValue is meant to be the previous value of counter.
580  //   corrected_diff = solve_ac_overflow(oldValueOfTheCounter, actualValue);  //
581  //   ...  // Example:
582  // }  // for(...) {
583  //  //   ...
584  //  //   corrected_diff = solve_ac_overflow(oldValueOfTheCounter, actualValue);
585  // Returns the corrected difference between counter and oldValue and  //   ...
586  // set oldValue to the value of counter.  // }
587  // Attention: oldValue is a reference.  //
588  Int_t solve_ac_overflow(Int_t& oldValue, Int_t counter)  //
589  {  // Returns the corrected difference between counter and oldValue and
590    Int_t truediff = 0;  // set oldValue to the value of counter.
591    // Attention: oldValue is a reference.
592    if (counter < oldValue)  // overflow!  Int_t solve_ac_overflow(Int_t& oldValue, Int_t counter)
593      truediff = 0xFFFF - oldValue + counter;  {
594    else    Int_t truediff = 0;
595      truediff = counter - oldValue;  
596      if (counter < oldValue)  // overflow!
597    oldValue = counter;      truediff = 0xFFFF - oldValue + counter;
598      else
599    return truediff;      truediff = counter - oldValue;
600  }  
601      oldValue = counter;
602    
603      return truediff;
604    }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.23