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

Diff of /quicklook/OrbitalRate/src/OrbitalRate.cpp

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

revision 1.1 by pam-rm2, Tue Dec 5 19:49:19 2006 UTC revision 1.9 by pam-rm2, Sun Nov 4 16:01:37 2007 UTC
# Line 50  int main(int argc, char* argv[]){ Line 50  int main(int argc, char* argv[]){
50    int offDate = 20060928;    int offDate = 20060928;
51    //  int offDate = 20060614;    //  int offDate = 20060614;
52    int offTime = 210000;    int offTime = 210000;
53      bool field = false;
54    
55    if (argc < 2){    if (argc < 2){
56      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 66  int main(int argc, char* argv[]){ Line 67  int main(int argc, char* argv[]){
67      printf( "\t -outDir[path]          Path where to put the output.\n");      printf( "\t -outDir[path]          Path where to put the output.\n");
68      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");
69      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");
70        printf( "\t -field               Produce maps of the magnetic field \n");
71      exit(1);      exit(1);
72    }    }
73    
# Line 79  int main(int argc, char* argv[]){ Line 81  int main(int argc, char* argv[]){
81    }    }
82    
83    for (int i = 2; i < argc; i++){    for (int i = 2; i < argc; i++){
84        if (!strcmp(argv[i], "-field")){
85          field = true;
86          i++;
87          continue;
88        }
89    
90      if (!strcmp(argv[i], "-outDir")){      if (!strcmp(argv[i], "-outDir")){
91        if (++i >= argc){        if (++i >= argc){
92          printf( "-outDir needs arguments. \n");          printf( "-outDir needs arguments. \n");
# Line 138  int main(int argc, char* argv[]){ Line 146  int main(int argc, char* argv[]){
146    }    }
147    
148    if (mapFile != ""){    if (mapFile != ""){
149      Rate(rootFile, outDir, mapFile, tleFile, offDate, offTime);      Rate(rootFile, outDir, mapFile, tleFile, offDate, offTime, field);
150    } else {    } else {
151      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");
152      printf("Try '--help' for more information. \n");      printf("Try '--help' for more information. \n");
# Line 207  void InitStyle() { Line 215  void InitStyle() {
215  }  }
216    
217    
218  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, bool field = false)
219  {  {
220    // **** Offset to temporarily correct the TDatime bug ****/    // **** Offset to temporarily correct the TDatime bug ****/
221    offTime += 10000;    //  offTime += 10000;
222    //********************************************************/    //********************************************************/
223    
224    TTree                  *tr         = 0;    TTree                  *tr         = 0;
# Line 272  void Rate(TString *filename, TString out Line 280  void Rate(TString *filename, TString out
280    // Magnetic field histograms.  I use always the suffix _counter    // Magnetic field histograms.  I use always the suffix _counter
281    // because they are not normalized.  Imagine that an instrument    // because they are not normalized.  Imagine that an instrument
282    // give us the value of the magnetic field for each event.    // give us the value of the magnetic field for each event.
283    TH2F *hbabs_counter = new TH2F("hbabs_counter", "B module", 360, -180, 180, 180, -90, 90);    TH2F *hbabs_counter;
284    TH2F *hbnorth_counter = new TH2F("hbnorth_counter", "B north", 360, -180, 180, 180, -90, 90);    TH2F *hbnorth_counter;
285    TH2F *hbdown_counter = new TH2F("hbdown_counter", "B down", 360, -180, 180, 180, -90, 90);    TH2F *hbdown_counter;
286    TH2F *hbeast_counter = new TH2F("hbeast_counter", "B east", 360, -180, 180, 180, -90, 90);    TH2F *hbeast_counter;
287    TH2F *hb0_counter = new TH2F("hb0_counter", "B_0", 360, -180, 180, 180, -90, 90);    TH2F *hb0_counter;
288    TH2F *hl_counter = new TH2F("hl_counter", "l", 360, -180, 180, 180, -90, 90);    TH2F *hl_counter;
289    
290      if(field) {
291        hbabs_counter = new TH2F("hbabs_counter", "B module", 360, -180, 180, 180, -90, 90);
292        hbnorth_counter = new TH2F("hbnorth_counter", "B north", 360, -180, 180, 180, -90, 90);
293        hbdown_counter = new TH2F("hbdown_counter", "B down", 360, -180, 180, 180, -90, 90);
294        hbeast_counter = new TH2F("hbeast_counter", "B east", 360, -180, 180, 180, -90, 90);
295        hb0_counter = new TH2F("hb0_counter", "B_0", 360, -180, 180, 180, -90, 90);
296        hl_counter = new TH2F("hl_counter", "l", 360, -180, 180, 180, -90, 90);
297      }
298    
299    // Get a char* to "file" from "/dir1/dir2/.../file.root"    // Get a char* to "file" from "/dir1/dir2/.../file.root"
300    TString basename;    TString basename;
# Line 307  void Rate(TString *filename, TString out Line 324  void Rate(TString *filename, TString out
324      exit(EXIT_FAILURE);      exit(EXIT_FAILURE);
325    }    }
326    
327    //Get the Julian date of the Resours offset    // Here I do: resurs offset + timesync
328    TDatime offRes = TDatime(offDate, offTime);    TDatime offRes = TDatime(offDate, offTime);
329    // Add to the Resours Offset the timesync.  This is now the date at    TTimeStamp offResTS = TTimeStamp(offRes.GetYear(), offRes.GetMonth(), offRes.GetDay(), offRes.GetHour(), offRes.GetMinute(), offRes.GetSecond(), 0, kTRUE, timesync);
330    // the moment of the timesync.  //   cout<<"timesync="<<timesync<<"  obt_timesync="<<obt_timesync<<endl;
331    offRes.Set(offRes.Convert() + (UInt_t) timesync);  //   cout<<"offResTS.GetSec()="<<offResTS.GetSec()<<endl;
332    //   cout<<"offResTS.AsString()="<<offResTS.AsString()<<endl<<endl;
333    
334    // Now I need a pointer to a cTle object.  The class misses a    // Now I need a pointer to a cTle object.  The class misses a
335    // constructor without arguments, so we have to give it a dummy TLE.    // constructor without arguments, so we have to give it a dummy TLE.
# Line 322  void Rate(TString *filename, TString out Line 340  void Rate(TString *filename, TString out
340    
341    // If we have to use a TLE file, call getTle().    // If we have to use a TLE file, call getTle().
342    if (tleFile != "")    if (tleFile != "")
343      tle1 = getTle(tleFile, offRes);      tle1 = getTle(tleFile, offResTS); // modify getTle() to use offResTS!
344      else
345        cout<<"OrbitalRate: Warning!!! No tle file supplied.\n";
346    
347      // Here I do: resurs offset + timesync - obt of the timesync
348      //  cout<<"here0 "<<offResTS.GetSec()<<endl;
349      //  offResTS.Set(offResTS.GetSec() - obt_timesync, kTRUE, 0, kFALSE);
350      offResTS.Set(offResTS.GetSec() - obt_timesync, kFALSE, 0, kFALSE);
351      //  cout<<"here1 "<<offResTS.GetSec()<<endl;
352    
353    cOrbit       orbit(*tle1);    cOrbit       orbit(*tle1);
354    cEci         eci;    cEci         eci;
355    cCoordGeo    coo;    cCoordGeo    coo;
356    
357    // offRes is now "offset date" + timesync.  Now I subtract the obt    // Here I do: resurs offset + timesync - obt of the timesync - tle time
358    // of the timesync.  Remember that the time of the event from the    TTimeStamp tledate = getTleDatetime(tle1);
   // tle date is:  
   // tle date - (offset date + timesync - obt timesync + obt event).  
   offRes.Set(offRes.Convert() - (UInt_t) obt_timesync);  
   
   // Get the Julian date of the TLE epoch  
   string datetime = getTleDatetime(tle1);  
   TDatime tledate = TDatime(datetime.c_str());  
   
359    cJulian jdatetime = cJulian((int) (tle1->getField(cTle::FLD_EPOCHYEAR)+2e3), tle1->getField(cTle::FLD_EPOCHDAY));    cJulian jdatetime = cJulian((int) (tle1->getField(cTle::FLD_EPOCHYEAR)+2e3), tle1->getField(cTle::FLD_EPOCHDAY));
360    int pYear, pMon; double pDOM;    int pYear, pMon; double pDOM;
361    jdatetime.getComponent(&pYear, &pMon, &pDOM);    jdatetime.getComponent(&pYear, &pMon, &pDOM);
362      offsetTime = ((Long64_t) offResTS.GetSec() - (Long64_t) tledate.GetSec());
363    offsetTime = ((Long64_t) offRes.Convert() - (Long64_t) tledate.Convert());  //   cerr<<"offsetTime="<<offsetTime<<endl
364    //       <<"offResTS.GetSec()="<<offResTS.GetSec()<<endl
365    //       <<"tledate.GetSec()="<<tledate.GetSec()<<endl<<endl;
366    
367    /********** Magnetic Field **************/    /********** Magnetic Field **************/
368    // Check that all this is correct!    // Check that all this is correct!
369    float dimo = 0.0; // dipole moment (computed from dat files)    float br, btheta, bphi;
   float bnorth, beast, bdown, babs;  
   float xl; // L value  
   float icode; // code value for L accuracy (see fortran code)  
   float bab1; // What's  the difference with babs?  
   float stps = 0.005; // step size for field line tracing  
   float bdel = 0.01; // required accuracy  
   float bequ;  // equatorial b value (also called b_0)  
   bool value = 0; // false if bequ is not the minimum b value  
   float rr0; // equatorial radius normalized to earth radius  
   
   // Initialize fortran routines!!!  
   initize_();  
370    
371    // I can now compute the magnetic dipole moment at the actual date,    // I can now compute the magnetic dipole moment at the actual date,
372    // using the cJulian date.  I don't to recompute it for every event    // using the cJulian date.  I don't to recompute it for every event
373    // beacause changes are not relevant at all.    // beacause changes are not relevant at all.
374    Int_t y = tledate.GetYear();    UInt_t y, m, d;
375    Int_t m = tledate.GetMonth();    tledate.GetDate(kTRUE, 0, &y, &m, &d);
   Int_t d = tledate.GetDay();  
376    float year = (float) y + (m*31+d)/365;    float year = (float) y + (m*31+d)/365;
377    
378    // Compute the magnetic dipole moment.    // Initialize common data for geopack
379    feldcof_(&year, &dimo);    if(field)
380        recalc_(y, m*31+d, 0, 0, 0);
381    /********** Magnetic Field **************/    /********** Magnetic Field **************/
382    
383    tr = (TTree*)rootFile->Get("Physics");    tr = (TTree*)rootFile->Get("Physics");
# Line 416  void Rate(TString *filename, TString out Line 424  void Rate(TString *filename, TString out
424        ph = eh->GetPscuHeader();        ph = eh->GetPscuHeader();
425    
426        // obt in ms        // obt in ms
427        ULong64_t obt = ph->GetOrbitalTime();        UInt_t obt = (UInt_t) ph->GetOrbitalTime();
428    
429        // timeElapsedFromTLE is the difference, in seconds, between the        // timeElapsedFromTLE is the difference, in seconds, between the
430        // event and the tle date.  I use seconds and not milliseconds        // event and the tle date.  I use seconds and not milliseconds
431        // because the indetermination on the timesync is about 1s.        // because the indetermination on the timesync is about 1s.
432        timeElapsedFromTLE = offsetTime + obt/1000;        timeElapsedFromTLE = offsetTime + obt/1000;
433          if(!i) cerr<<"1st event: timeElapsedFromTLE="<<timeElapsedFromTLE<<endl;
434    
435        // I also need the abstime in seconds rounded to the lower        // I also need the abstime in seconds rounded to the lower
436        // value.  Every second, we set a_second_is_over to true.  Only        // value.  Every second, we set a_second_is_over to true.  Only
# Line 469  void Rate(TString *filename, TString out Line 478  void Rate(TString *filename, TString out
478        alt = coo.m_Alt;        alt = coo.m_Alt;
479    
480        /********** Magnetic Field **************/        /********** Magnetic Field **************/
481        feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs);        if(field)
482        shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1);          igrf_geo__((coo.m_Alt+6371.2)/6371.2, M_PI/2.-coo.m_Lat, coo.m_Lon, br, btheta, bphi);
483        findb0_(&stps, &bdel, &value, &bequ, &rr0);        //      cout<<"("<<(coo.m_Alt+6371.2)/6371.2<<", "<<M_PI/2.-coo.m_Lat<<", "<<coo.m_Lon<<")"<<endl;
484        /********** Magnetic Field **************/        /********** Magnetic Field **************/
485    
486        // serve fare il controllo deltatime < 1?        // serve fare il controllo deltatime < 1?
487        if (deltaTime > 1) cout << endl << "******** deltaTime<1 ********" << endl;        if (deltaTime > 1) cout << endl << "******** deltaTime<1 ********" << endl;
488        // Does nothing for the first two events or if acquisition time if more        // Does nothing for the first two events or if acquisition time if more
489        // than 1s.        // than 1s.
490        if(!i || (deltaTime > 1)) continue;        if(i<1 || (deltaTime > 1)) continue;
491    
492        // CAS3 and CAS4 are not rates but only counters.  So I fill        // CAS3 and CAS4 are not rates but only counters.  So I fill
493        // with the bin with the difference beetween the actual counter        // with the bin with the difference beetween the actual counter
# Line 498  void Rate(TString *filename, TString out Line 507  void Rate(TString *filename, TString out
507        // this values but I need to count how many times I fill        // this values but I need to count how many times I fill
508        // each bin.  This is done by the histogram event_counter.        // each bin.  This is done by the histogram event_counter.
509        // I will normalize later.        // I will normalize later.
510        hbabs_counter->Fill(lon, lat, babs);        if(field) {
511        hbnorth_counter->Fill(lon, lat, bnorth);          hbabs_counter->Fill(lon, lat, sqrt(br*br+btheta*btheta+bphi*bphi)*1e-5);
512        hbdown_counter->Fill(lon, lat, bdown);          hbnorth_counter->Fill(lon, lat, -btheta*1e-5);
513        hbeast_counter->Fill(lon, lat, beast);          hbdown_counter->Fill(lon, lat, -br*1e-5);
514        hb0_counter->Fill(lon, lat, bequ);          hbeast_counter->Fill(lon, lat, bphi*1e-5);
515        hl_counter->Fill(lon, lat, xl);        }
     
516        // This histograms is now filled with the number of entries.        // This histograms is now filled with the number of entries.
517        // Below we will divide with the time (in seconds) to get        // Below we will divide with the time (in seconds) to get
518        // event rate per bin.        // event rate per bin.
# Line 570  void Rate(TString *filename, TString out Line 578  void Rate(TString *filename, TString out
578    TH2F *trigAndOr_rate = (TH2F*) trigAndOr_counter->Clone("trigAndOr_rate");    TH2F *trigAndOr_rate = (TH2F*) trigAndOr_counter->Clone("trigAndOr_rate");
579    TH2F *trigAndAnd_rate = (TH2F*) trigAndAnd_counter->Clone("trigAndAnd_rate");    TH2F *trigAndAnd_rate = (TH2F*) trigAndAnd_counter->Clone("trigAndAnd_rate");
580    TH2F *nd_rate = (TH2F*) nd_counter->Clone("nd_rate");    TH2F *nd_rate = (TH2F*) nd_counter->Clone("nd_rate");
581    TH2F *hbabs_norm = (TH2F*) hbabs_counter->Clone("hbabs_norm");  
582    TH2F *hbnorth_norm = (TH2F*) hbnorth_counter->Clone("hbnorth_norm");    TH2F *hbabs_norm;
583    TH2F *hbdown_norm = (TH2F*) hbabs_counter->Clone("hbdown_norm");    TH2F *hbnorth_norm;
584    TH2F *hbeast_norm = (TH2F*) hbabs_counter->Clone("hbeast_norm");    TH2F *hbdown_norm;
585    TH2F *hb0_norm = (TH2F*) hb0_counter->Clone("hb0_norm");    TH2F *hbeast_norm;
586    TH2F *hl_norm = (TH2F*) hl_counter->Clone("hl_norm");  
587      if(field) {
588        hbabs_norm = (TH2F*) hbabs_counter->Clone("hbabs_norm");
589        hbnorth_norm = (TH2F*) hbnorth_counter->Clone("hbnorth_norm");
590        hbdown_norm = (TH2F*) hbabs_counter->Clone("hbdown_norm");
591        hbeast_norm = (TH2F*) hbabs_counter->Clone("hbeast_norm");
592      }
593    
594    // Now we divide each histogram _counter with the time histogram    // Now we divide each histogram _counter with the time histogram
595    // obtBinTime to have an histogram _rate.  Note that, when a second    // obtBinTime to have an histogram _rate.  Note that, when a second
# Line 591  void Rate(TString *filename, TString out Line 605  void Rate(TString *filename, TString out
605    trigS111A_rate->Divide(trigS111A_counter, obtBinTime, 1, 1, "");    trigS111A_rate->Divide(trigS111A_counter, obtBinTime, 1, 1, "");
606    oss.str("");    oss.str("");
607    oss << basename.Data() << "_orbit_trigS111A.png";    oss << basename.Data() << "_orbit_trigS111A.png";
608      trigS111A_rate->SetMinimum(9);
609    printHist(trigS111A_rate, mapFile, outDirectory, oss.str().c_str(), "S111A (Hz)", -width, height, true, 0);    printHist(trigS111A_rate, mapFile, outDirectory, oss.str().c_str(), "S111A (Hz)", -width, height, true, 0);
610    
611    antiCAS4_rate->Divide(antiCAS4_counter, obtBinTime, 1, 1, "");    antiCAS4_rate->Divide(antiCAS4_counter, obtBinTime, 1, 1, "");
612    oss.str("");    oss.str("");
613    oss << basename.Data() << "_orbit_CAS4.png";    oss << basename.Data() << "_orbit_CAS4.png";
614      antiCAS4_rate->SetMinimum(99);
615    printHist(antiCAS4_rate, mapFile, outDirectory, oss.str().c_str(), "CAS4 (Hz)", -width, height, true, 0);    printHist(antiCAS4_rate, mapFile, outDirectory, oss.str().c_str(), "CAS4 (Hz)", -width, height, true, 0);
616    
617    antiCAS3_rate->Divide(antiCAS3_counter, obtBinTime, 1, 1, "");    antiCAS3_rate->Divide(antiCAS3_counter, obtBinTime, 1, 1, "");
618    oss.str("");    oss.str("");
619    oss << basename.Data() << "_orbit_CAS3.png";    oss << basename.Data() << "_orbit_CAS3.png";
620      antiCAS3_rate->SetMinimum(99);
621    printHist(antiCAS3_rate, mapFile, outDirectory, oss.str().c_str(), "CAS3 (Hz)", -width, height, true, 0);    printHist(antiCAS3_rate, mapFile, outDirectory, oss.str().c_str(), "CAS3 (Hz)", -width, height, true, 0);
622    
623    event_rate->Divide(event_counter, obtBinTime, 1, 1, "");    event_rate->Divide(event_counter, obtBinTime, 1, 1, "");
# Line 611  void Rate(TString *filename, TString out Line 628  void Rate(TString *filename, TString out
628    trigS11andS12_rate->Divide(trigS11andS12_counter, obtBinTime, 1, 1, "");    trigS11andS12_rate->Divide(trigS11andS12_counter, obtBinTime, 1, 1, "");
629    oss.str("");    oss.str("");
630    oss << basename.Data() << "_orbit_trigS11andS12.png";    oss << basename.Data() << "_orbit_trigS11andS12.png";
631      trigS11andS12_rate->SetMinimum(99);
632    printHist(trigS11andS12_rate, mapFile, outDirectory, oss.str().c_str(), "(S11*S12) (Hz)", -width, height, 1, 0);    printHist(trigS11andS12_rate, mapFile, outDirectory, oss.str().c_str(), "(S11*S12) (Hz)", -width, height, 1, 0);
633    
634    trigS12andS21andS22_rate->Divide(trigS12andS21andS22_counter, obtBinTime, 1, 1, "");    trigS12andS21andS22_rate->Divide(trigS12andS21andS22_counter, obtBinTime, 1, 1, "");
635    oss.str("");    oss.str("");
636    oss << basename.Data() << "_orbit_trigS12andS21andS22.png";    oss << basename.Data() << "_orbit_trigS12andS21andS22.png";
637      trigS12andS21andS22_rate->SetMinimum(9);
638    printHist(trigS12andS21andS22_rate, mapFile, outDirectory, oss.str().c_str(), "(S12*S12*S21) (Hz)", -width, height, true, 0);    printHist(trigS12andS21andS22_rate, mapFile, outDirectory, oss.str().c_str(), "(S12*S12*S21) (Hz)", -width, height, true, 0);
639    
640    trigAndOr_rate->Divide(trigAndOr_counter, obtBinTime, 1, 1, "");    trigAndOr_rate->Divide(trigAndOr_counter, obtBinTime, 1, 1, "");
# Line 637  void Rate(TString *filename, TString out Line 656  void Rate(TString *filename, TString out
656    // fill the bins with the values of the magnetic field for each    // fill the bins with the values of the magnetic field for each
657    // event, we need to divide with the number of fills done, that is    // event, we need to divide with the number of fills done, that is
658    // event_counter.    // event_counter.
659    hbabs_norm->Divide(hbabs_counter, event_counter, 1, 1, "");    if(field) {
660    oss.str("");      hbabs_norm->Divide(hbabs_counter, event_counter, 1, 1, "");
661    oss << basename.Data() << "_orbit_Babs.png";      oss.str("");
662    printHist(hbabs_norm, mapFile, outDirectory, oss.str().c_str(), "B abs (G)", -width, height, 0, 0);      oss << basename.Data() << "_orbit_Babs.png";
663        printHist(hbabs_norm, mapFile, outDirectory, oss.str().c_str(), "B abs (G)", -width, height, 0, 0);
664    hbnorth_norm->Divide(hbnorth_counter, event_counter, 1, 1, "");  
665    oss.str("");      hbnorth_norm->Divide(hbnorth_counter, event_counter, 1, 1, "");
666    oss << basename.Data() << "_orbit_Bnorth.png";      oss.str("");
667    printHist(hbnorth_norm, mapFile, outDirectory, oss.str().c_str(), "B north (G)", -width, height, 0, 1);      oss << basename.Data() << "_orbit_Bnorth.png";
668        printHist(hbnorth_norm, mapFile, outDirectory, oss.str().c_str(), "B north (G)", -width, height, 0, 1);
669    hbdown_norm->Divide(hbdown_counter, event_counter, 1, 1, "");  
670    oss.str("");      hbdown_norm->Divide(hbdown_counter, event_counter, 1, 1, "");
671    oss << basename.Data() << "_orbit_Bdown.png";      oss.str("");
672    printHist(hbdown_norm, mapFile, outDirectory, oss.str().c_str(), "B down (G)", -width, height, 0, 1);      oss << basename.Data() << "_orbit_Bdown.png";
673        printHist(hbdown_norm, mapFile, outDirectory, oss.str().c_str(), "B down (G)", -width, height, 0, 1);
674    hbeast_norm->Divide(hbeast_counter, event_counter, 1, 1, "");  
675    oss.str("");      hbeast_norm->Divide(hbeast_counter, event_counter, 1, 1, "");
676    oss << basename.Data() << "_orbit_Beast.png";      oss.str("");
677    printHist(hbeast_norm, mapFile, outDirectory, oss.str().c_str(), "B east (G)", -width, height, 0, 1);      oss << basename.Data() << "_orbit_Beast.png";
678        printHist(hbeast_norm, mapFile, outDirectory, oss.str().c_str(), "B east (G)", -width, height, 0, 1);
679    hb0_norm->Divide(hb0_counter, event_counter, 1, 1, "");    }
   oss.str("");  
   oss << basename.Data() << "_orbit_B0.png";  
   printHist(hb0_norm, mapFile, outDirectory, oss.str().c_str(), "B_0 (G)", -width, height, 0, 0);  
   
   hl_norm->Divide(hl_counter, event_counter, 1, 1, "");  
   oss.str("");  
   oss << basename.Data() << "_orbit_L.png";  
   printHist(hl_norm, mapFile, outDirectory, oss.str().c_str(), "L shell", -width, height, 0, 0);  
   
680    
681    delete obtBinTime;    delete obtBinTime;
682    delete event_counter;    delete event_counter;
# Line 690  void Rate(TString *filename, TString out Line 700  void Rate(TString *filename, TString out
700    delete trigS111A_rate;    delete trigS111A_rate;
701    delete trigS12andS21andS22_rate;    delete trigS12andS21andS22_rate;
702    
703    delete hbabs_counter;    if(field) {
704    delete hbnorth_counter;      delete hbabs_counter;
705    delete hbdown_counter;      delete hbnorth_counter;
706    delete hbeast_counter;      delete hbdown_counter;
707    delete hb0_counter;      delete hbeast_counter;
708    delete hl_counter;      delete hbabs_norm;
709    delete hbabs_norm;      delete hbnorth_norm;
710    delete hbnorth_norm;      delete hbdown_norm;
711    delete hbdown_norm;      delete hbeast_norm;
712    delete hbeast_norm;    }
   delete hb0_norm;  
   delete hl_norm;  
713    
714    rootFile->Close();    rootFile->Close();
715  }  }
# Line 725  void Rate(TString *filename, TString out Line 733  void Rate(TString *filename, TString out
733  // Scale() and Merge()).  // Scale() and Merge()).
734  //  //
735  // This function depends on InitStyle();  // This function depends on InitStyle();
736  int printHist(TH2F *h, TString mapFile, TString outDirectory, TString outputFilename, char *title, int width, int height, bool use_log, bool bool_shift)  int printHist(TH2F *h, TString mapFile, TString outDirectory, TString outputFilename, const char *title, int width, int height, bool use_log, bool bool_shift)
737  {  {
738    InitStyle();    InitStyle();
739    
740    // Create a canvas and draw the TH2F with a nice colormap for z    // Create a canvas and draw the TH2F with a nice colormap for z
741    // values, using log scale for z values, if requested, and setting    // values, using log scale for z values, if requested, and setting
742    // some title.    // some title.
743    TCanvas *canvas = new TCanvas("h", "passed histogram", width*2, height*2);    TCanvas *canvas = new TCanvas("h", "h histogram", width*2, height*2);
744    
745    if(use_log) {    if(use_log) canvas->SetLogz();
     h->SetMinimum(1);  
     canvas->SetLogz();  
   }  
746    
747    h->SetTitle(title);    h->SetTitle(title);
748    h->SetXTitle("Longitude (deg)");    h->SetXTitle("Longitude (deg)");
# Line 826  void saveHist(TH1 *h, TString savetoroot Line 831  void saveHist(TH1 *h, TString savetoroot
831  // querying the database with the RESURS DK-1 id number 29228,  // querying the database with the RESURS DK-1 id number 29228,
832  // selecting the widest timespan, including the satellite name in the  // selecting the widest timespan, including the satellite name in the
833  // results.  // results.
834  cTle *getTle(TString tleFile, TDatime offRes)  cTle *getTle(TString tleFile, TTimeStamp offResTS)
835  {  {
836    Float_t tledatefromfile, tledatefromroot;    Float_t tledatefromfile, tledatefromroot;
837    fstream tlefile(tleFile.Data(), ios::in);    fstream tlefile(tleFile.Data(), ios::in);
# Line 856  cTle *getTle(TString tleFile, TDatime of Line 861  cTle *getTle(TString tleFile, TDatime of
861    // Sort by date    // Sort by date
862    sort(ctles.begin(), ctles.end(), compTLE);    sort(ctles.begin(), ctles.end(), compTLE);
863    
864    tledatefromroot = (offRes.GetYear()-2000)*1e3 + (offRes.Convert() - (TDatime(offRes.GetYear(), 1, 1, 0, 0, 0)).Convert())/ (24.*3600.);    UInt_t year, month, day;
865      offResTS.GetDate(kTRUE, 0, &year, &month, &day);
866      TTimeStamp firstofjan = TTimeStamp(year, 1, 1, 0, 0, 0);
867      tledatefromroot = (year-2000)*1e3 + (offResTS.GetSec() - firstofjan.GetSec())/(24.*3600.);
868    
869    for(iter = ctles.begin(); iter != ctles.end(); iter++) {    for(iter = ctles.begin(); iter != ctles.end(); iter++) {
870      cTle *tle = *iter;      cTle *tle = *iter;
# Line 900  float getTleJulian(cTle *tle) { Line 908  float getTleJulian(cTle *tle) {
908    
909  // Look for a timesync in the TFile rootFile.  Set timesync and  // Look for a timesync in the TFile rootFile.  Set timesync and
910  // obt_timesync.  Returns 1 if timesync is found, 0 otherwise.  // obt_timesync.  Returns 1 if timesync is found, 0 otherwise.
911  int lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) {  UInt_t lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) {
912    *timesync = -1;  // will be != -1 if found    *timesync = -1;  // will be != -1 if found
913    
914    ULong64_t             nevents    = 0;    ULong64_t             nevents    = 0;
# Line 971  TH2F* shiftHist(TH2F* h, Float_t shift) Line 979  TH2F* shiftHist(TH2F* h, Float_t shift)
979  }  }
980    
981    
982  // Return a string like YYYY-MM-DD hh:mm:ss, a datetime format.  //
983  string getTleDatetime(cTle *tle)  // Returns the tle date as a TTimeStamp object.
984    //
985    TTimeStamp getTleDatetime(cTle *tle)
986  {  {
987    int year, mon, day, hh, mm, ss;    int year, mon, day, hh, mm, ss;
988    double dom; // day of month (is double!)    double dom; // day of month (is double!)
# Line 992  string getTleDatetime(cTle *tle) Line 1002  string getTleDatetime(cTle *tle)
1002    ss = (int) floor( ((((dom - day) * 24 - hh) * 60 - mm) * 60));    ss = (int) floor( ((((dom - day) * 24 - hh) * 60 - mm) * 60));
1003    //  ms = (int) floor( (((((dom - day) * 24 - hh) * 60 - mm) * 60) - ss) * 1000);    //  ms = (int) floor( (((((dom - day) * 24 - hh) * 60 - mm) * 60) - ss) * 1000);
1004    
1005    date << year << "-" << mon << "-" << day << " " << hh << ":" << mm << ":" << ss;    TTimeStamp t = TTimeStamp(year, mon, day, hh, mm, ss, 0, true);
1006    
1007    return date.str();    return t;
1008  }  }
1009    
1010  //  //
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.23