/[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.4 by pam-rm2, Thu Dec 7 13:27:30 2006 UTC revision 1.7 by pam-rm2, Mon Apr 2 19:44:38 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 357  void Rate(TString *filename, TString out Line 374  void Rate(TString *filename, TString out
374    float year = (float) y + (m*31+d)/365;    float year = (float) y + (m*31+d)/365;
375    
376    // Initialize common data for geopack    // Initialize common data for geopack
377    recalc_(y, m*31+d, 0, 0, 0);    if(field)
378        recalc_(y, m*31+d, 0, 0, 0);
379    /********** Magnetic Field **************/    /********** Magnetic Field **************/
380    
381    tr = (TTree*)rootFile->Get("Physics");    tr = (TTree*)rootFile->Get("Physics");
# Line 404  void Rate(TString *filename, TString out Line 422  void Rate(TString *filename, TString out
422        ph = eh->GetPscuHeader();        ph = eh->GetPscuHeader();
423    
424        // obt in ms        // obt in ms
425        ULong64_t obt = ph->GetOrbitalTime();        UInt_t obt = (UInt_t) ph->GetOrbitalTime();
426    
427        // timeElapsedFromTLE is the difference, in seconds, between the        // timeElapsedFromTLE is the difference, in seconds, between the
428        // event and the tle date.  I use seconds and not milliseconds        // event and the tle date.  I use seconds and not milliseconds
# Line 457  void Rate(TString *filename, TString out Line 475  void Rate(TString *filename, TString out
475        alt = coo.m_Alt;        alt = coo.m_Alt;
476    
477        /********** Magnetic Field **************/        /********** Magnetic Field **************/
478        igrf_geo__((coo.m_Alt+6371.2)/6371.2, M_PI/2.-coo.m_Lat, coo.m_Lon, br, btheta, bphi);        if(field)
479            igrf_geo__((coo.m_Alt+6371.2)/6371.2, M_PI/2.-coo.m_Lat, coo.m_Lon, br, btheta, bphi);
480        //      cout<<"("<<(coo.m_Alt+6371.2)/6371.2<<", "<<M_PI/2.-coo.m_Lat<<", "<<coo.m_Lon<<")"<<endl;        //      cout<<"("<<(coo.m_Alt+6371.2)/6371.2<<", "<<M_PI/2.-coo.m_Lat<<", "<<coo.m_Lon<<")"<<endl;
481        /********** Magnetic Field **************/        /********** Magnetic Field **************/
482    
# Line 485  void Rate(TString *filename, TString out Line 504  void Rate(TString *filename, TString out
504        // this values but I need to count how many times I fill        // this values but I need to count how many times I fill
505        // each bin.  This is done by the histogram event_counter.        // each bin.  This is done by the histogram event_counter.
506        // I will normalize later.        // I will normalize later.
507        hbabs_counter->Fill(lon, lat, sqrt(br*br+btheta*btheta+bphi*bphi)*1e-5);        if(field) {
508        hbnorth_counter->Fill(lon, lat, -btheta*1e-5);          hbabs_counter->Fill(lon, lat, sqrt(br*br+btheta*btheta+bphi*bphi)*1e-5);
509        hbdown_counter->Fill(lon, lat, -br*1e-5);          hbnorth_counter->Fill(lon, lat, -btheta*1e-5);
510        hbeast_counter->Fill(lon, lat, bphi*1e-5);          hbdown_counter->Fill(lon, lat, -br*1e-5);
511              hbeast_counter->Fill(lon, lat, bphi*1e-5);
512          }
513        // This histograms is now filled with the number of entries.        // This histograms is now filled with the number of entries.
514        // Below we will divide with the time (in seconds) to get        // Below we will divide with the time (in seconds) to get
515        // event rate per bin.        // event rate per bin.
# Line 555  void Rate(TString *filename, TString out Line 575  void Rate(TString *filename, TString out
575    TH2F *trigAndOr_rate = (TH2F*) trigAndOr_counter->Clone("trigAndOr_rate");    TH2F *trigAndOr_rate = (TH2F*) trigAndOr_counter->Clone("trigAndOr_rate");
576    TH2F *trigAndAnd_rate = (TH2F*) trigAndAnd_counter->Clone("trigAndAnd_rate");    TH2F *trigAndAnd_rate = (TH2F*) trigAndAnd_counter->Clone("trigAndAnd_rate");
577    TH2F *nd_rate = (TH2F*) nd_counter->Clone("nd_rate");    TH2F *nd_rate = (TH2F*) nd_counter->Clone("nd_rate");
578    TH2F *hbabs_norm = (TH2F*) hbabs_counter->Clone("hbabs_norm");  
579    TH2F *hbnorth_norm = (TH2F*) hbnorth_counter->Clone("hbnorth_norm");    TH2F *hbabs_norm;
580    TH2F *hbdown_norm = (TH2F*) hbabs_counter->Clone("hbdown_norm");    TH2F *hbnorth_norm;
581    TH2F *hbeast_norm = (TH2F*) hbabs_counter->Clone("hbeast_norm");    TH2F *hbdown_norm;
582      TH2F *hbeast_norm;
583    
584      if(field) {
585        hbabs_norm = (TH2F*) hbabs_counter->Clone("hbabs_norm");
586        hbnorth_norm = (TH2F*) hbnorth_counter->Clone("hbnorth_norm");
587        hbdown_norm = (TH2F*) hbabs_counter->Clone("hbdown_norm");
588        hbeast_norm = (TH2F*) hbabs_counter->Clone("hbeast_norm");
589      }
590    
591    // Now we divide each histogram _counter with the time histogram    // Now we divide each histogram _counter with the time histogram
592    // obtBinTime to have an histogram _rate.  Note that, when a second    // obtBinTime to have an histogram _rate.  Note that, when a second
# Line 625  void Rate(TString *filename, TString out Line 653  void Rate(TString *filename, TString out
653    // fill the bins with the values of the magnetic field for each    // fill the bins with the values of the magnetic field for each
654    // 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
655    // event_counter.    // event_counter.
656    hbabs_norm->Divide(hbabs_counter, event_counter, 1, 1, "");    if(field) {
657    oss.str("");      hbabs_norm->Divide(hbabs_counter, event_counter, 1, 1, "");
658    oss << basename.Data() << "_orbit_Babs.png";      oss.str("");
659    printHist(hbabs_norm, mapFile, outDirectory, oss.str().c_str(), "B abs (G)", -width, height, 0, 0);      oss << basename.Data() << "_orbit_Babs.png";
660        printHist(hbabs_norm, mapFile, outDirectory, oss.str().c_str(), "B abs (G)", -width, height, 0, 0);
661    hbnorth_norm->Divide(hbnorth_counter, event_counter, 1, 1, "");  
662    oss.str("");      hbnorth_norm->Divide(hbnorth_counter, event_counter, 1, 1, "");
663    oss << basename.Data() << "_orbit_Bnorth.png";      oss.str("");
664    printHist(hbnorth_norm, mapFile, outDirectory, oss.str().c_str(), "B north (G)", -width, height, 0, 1);      oss << basename.Data() << "_orbit_Bnorth.png";
665        printHist(hbnorth_norm, mapFile, outDirectory, oss.str().c_str(), "B north (G)", -width, height, 0, 1);
666    hbdown_norm->Divide(hbdown_counter, event_counter, 1, 1, "");  
667    oss.str("");      hbdown_norm->Divide(hbdown_counter, event_counter, 1, 1, "");
668    oss << basename.Data() << "_orbit_Bdown.png";      oss.str("");
669    printHist(hbdown_norm, mapFile, outDirectory, oss.str().c_str(), "B down (G)", -width, height, 0, 1);      oss << basename.Data() << "_orbit_Bdown.png";
670        printHist(hbdown_norm, mapFile, outDirectory, oss.str().c_str(), "B down (G)", -width, height, 0, 1);
671    hbeast_norm->Divide(hbeast_counter, event_counter, 1, 1, "");  
672    oss.str("");      hbeast_norm->Divide(hbeast_counter, event_counter, 1, 1, "");
673    oss << basename.Data() << "_orbit_Beast.png";      oss.str("");
674    printHist(hbeast_norm, mapFile, outDirectory, oss.str().c_str(), "B east (G)", -width, height, 0, 1);      oss << basename.Data() << "_orbit_Beast.png";
675        printHist(hbeast_norm, mapFile, outDirectory, oss.str().c_str(), "B east (G)", -width, height, 0, 1);
676      }
677    
678    delete obtBinTime;    delete obtBinTime;
679    delete event_counter;    delete event_counter;
# Line 668  void Rate(TString *filename, TString out Line 697  void Rate(TString *filename, TString out
697    delete trigS111A_rate;    delete trigS111A_rate;
698    delete trigS12andS21andS22_rate;    delete trigS12andS21andS22_rate;
699    
700    delete hbabs_counter;    if(field) {
701    delete hbnorth_counter;      delete hbabs_counter;
702    delete hbdown_counter;      delete hbnorth_counter;
703    delete hbeast_counter;      delete hbdown_counter;
704    delete hbabs_norm;      delete hbeast_counter;
705    delete hbnorth_norm;      delete hbabs_norm;
706    delete hbdown_norm;      delete hbnorth_norm;
707    delete hbeast_norm;      delete hbdown_norm;
708        delete hbeast_norm;
709      }
710    
711    rootFile->Close();    rootFile->Close();
712  }  }
# Line 871  float getTleJulian(cTle *tle) { Line 902  float getTleJulian(cTle *tle) {
902    
903  // Look for a timesync in the TFile rootFile.  Set timesync and  // Look for a timesync in the TFile rootFile.  Set timesync and
904  // obt_timesync.  Returns 1 if timesync is found, 0 otherwise.  // obt_timesync.  Returns 1 if timesync is found, 0 otherwise.
905  int lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) {  UInt_t lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) {
906    *timesync = -1;  // will be != -1 if found    *timesync = -1;  // will be != -1 if found
907    
908    ULong64_t             nevents    = 0;    ULong64_t             nevents    = 0;
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23