--- quicklook/OrbitalRate/src/OrbitalRate.cpp 2006/12/07 13:27:30 1.4 +++ quicklook/OrbitalRate/src/OrbitalRate.cpp 2007/11/03 22:11:18 1.8 @@ -50,6 +50,7 @@ int offDate = 20060928; // int offDate = 20060614; int offTime = 210000; + bool field = false; if (argc < 2){ printf("You have to insert at least the file to analyze and the mapFile \n"); @@ -66,6 +67,7 @@ printf( "\t -outDir[path] Path where to put the output.\n"); printf( "\t -offDate Date of resetting of the Resource counter [format YYMMDD (UTC date) default 20060928] \n"); printf( "\t -offTime Time of resetting of the Resource counter [format HHMMSS (UTC date) default 210000] \n"); + printf( "\t -field Produce maps of the magnetic field \n"); exit(1); } @@ -79,6 +81,12 @@ } for (int i = 2; i < argc; i++){ + if (!strcmp(argv[i], "-field")){ + field = true; + i++; + continue; + } + if (!strcmp(argv[i], "-outDir")){ if (++i >= argc){ printf( "-outDir needs arguments. \n"); @@ -138,7 +146,7 @@ } if (mapFile != ""){ - Rate(rootFile, outDir, mapFile, tleFile, offDate, offTime); + Rate(rootFile, outDir, mapFile, tleFile, offDate, offTime, field); } else { printf("You have to insert at least the file to analyze and the mapFile \n"); printf("Try '--help' for more information. \n"); @@ -207,10 +215,10 @@ } -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) { // **** Offset to temporarily correct the TDatime bug ****/ - offTime += 10000; + // offTime += 10000; //********************************************************/ TTree *tr = 0; @@ -272,12 +280,21 @@ // Magnetic field histograms. I use always the suffix _counter // because they are not normalized. Imagine that an instrument // give us the value of the magnetic field for each event. - TH2F *hbabs_counter = new TH2F("hbabs_counter", "B module", 360, -180, 180, 180, -90, 90); - TH2F *hbnorth_counter = new TH2F("hbnorth_counter", "B north", 360, -180, 180, 180, -90, 90); - TH2F *hbdown_counter = new TH2F("hbdown_counter", "B down", 360, -180, 180, 180, -90, 90); - TH2F *hbeast_counter = new TH2F("hbeast_counter", "B east", 360, -180, 180, 180, -90, 90); - TH2F *hb0_counter = new TH2F("hb0_counter", "B_0", 360, -180, 180, 180, -90, 90); - TH2F *hl_counter = new TH2F("hl_counter", "l", 360, -180, 180, 180, -90, 90); + TH2F *hbabs_counter; + TH2F *hbnorth_counter; + TH2F *hbdown_counter; + TH2F *hbeast_counter; + TH2F *hb0_counter; + TH2F *hl_counter; + + if(field) { + hbabs_counter = new TH2F("hbabs_counter", "B module", 360, -180, 180, 180, -90, 90); + hbnorth_counter = new TH2F("hbnorth_counter", "B north", 360, -180, 180, 180, -90, 90); + hbdown_counter = new TH2F("hbdown_counter", "B down", 360, -180, 180, 180, -90, 90); + hbeast_counter = new TH2F("hbeast_counter", "B east", 360, -180, 180, 180, -90, 90); + hb0_counter = new TH2F("hb0_counter", "B_0", 360, -180, 180, 180, -90, 90); + hl_counter = new TH2F("hl_counter", "l", 360, -180, 180, 180, -90, 90); + } // Get a char* to "file" from "/dir1/dir2/.../file.root" TString basename; @@ -307,11 +324,9 @@ exit(EXIT_FAILURE); } - //Get the Julian date of the Resours offset + // Here I do: resurs offset + timesync TDatime offRes = TDatime(offDate, offTime); - // Add to the Resours Offset the timesync. This is now the date at - // the moment of the timesync. - offRes.Set(offRes.Convert() + (UInt_t) timesync); + TTimeStamp offResTS = TTimeStamp(offRes.GetYear(), offRes.GetMonth(), offRes.GetDay(), offRes.GetHour(), offRes.GetMinute(), offRes.GetSecond(), 0, kTRUE, timesync); // Now I need a pointer to a cTle object. The class misses a // constructor without arguments, so we have to give it a dummy TLE. @@ -322,27 +337,23 @@ // If we have to use a TLE file, call getTle(). if (tleFile != "") - tle1 = getTle(tleFile, offRes); + tle1 = getTle(tleFile, offResTS); // modify getTle() to use offResTS! + else + cout<<"OrbitalRate: Warning!!! No tle file supplied.\n"; + + // Here I do: resurs offset + timesync - obt of the timesync + offResTS.Set(offResTS.GetSec() - obt_timesync, kTRUE, 0, kFALSE); cOrbit orbit(*tle1); cEci eci; cCoordGeo coo; - // offRes is now "offset date" + timesync. Now I subtract the obt - // of the timesync. Remember that the time of the event from the - // 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()); - + // Here I do: resurs offset + timesync - obt of the timesync - tle time + TTimeStamp tledate = getTleDatetime(tle1); cJulian jdatetime = cJulian((int) (tle1->getField(cTle::FLD_EPOCHYEAR)+2e3), tle1->getField(cTle::FLD_EPOCHDAY)); int pYear, pMon; double pDOM; jdatetime.getComponent(&pYear, &pMon, &pDOM); - - offsetTime = ((Long64_t) offRes.Convert() - (Long64_t) tledate.Convert()); + offsetTime = ((Long64_t) offResTS.GetSec() - (Long64_t) tledate.GetSec()); /********** Magnetic Field **************/ // Check that all this is correct! @@ -351,13 +362,16 @@ // I can now compute the magnetic dipole moment at the actual date, // using the cJulian date. I don't to recompute it for every event // beacause changes are not relevant at all. - Int_t y = tledate.GetYear(); - Int_t m = tledate.GetMonth(); - Int_t d = tledate.GetDay(); +// Int_t y = tledate.GetYear(); +// Int_t m = tledate.GetMonth(); +// Int_t d = tledate.GetDay(); + UInt_t y, m, d; + tledate.GetDate(kTRUE, 0, &y, &m, &d); float year = (float) y + (m*31+d)/365; // Initialize common data for geopack - recalc_(y, m*31+d, 0, 0, 0); + if(field) + recalc_(y, m*31+d, 0, 0, 0); /********** Magnetic Field **************/ tr = (TTree*)rootFile->Get("Physics"); @@ -404,7 +418,7 @@ ph = eh->GetPscuHeader(); // obt in ms - ULong64_t obt = ph->GetOrbitalTime(); + UInt_t obt = (UInt_t) ph->GetOrbitalTime(); // timeElapsedFromTLE is the difference, in seconds, between the // event and the tle date. I use seconds and not milliseconds @@ -457,7 +471,8 @@ alt = coo.m_Alt; /********** Magnetic Field **************/ - igrf_geo__((coo.m_Alt+6371.2)/6371.2, M_PI/2.-coo.m_Lat, coo.m_Lon, br, btheta, bphi); + if(field) + igrf_geo__((coo.m_Alt+6371.2)/6371.2, M_PI/2.-coo.m_Lat, coo.m_Lon, br, btheta, bphi); // cout<<"("<<(coo.m_Alt+6371.2)/6371.2<<", "<Fill(lon, lat, sqrt(br*br+btheta*btheta+bphi*bphi)*1e-5); - hbnorth_counter->Fill(lon, lat, -btheta*1e-5); - hbdown_counter->Fill(lon, lat, -br*1e-5); - hbeast_counter->Fill(lon, lat, bphi*1e-5); - + if(field) { + hbabs_counter->Fill(lon, lat, sqrt(br*br+btheta*btheta+bphi*bphi)*1e-5); + hbnorth_counter->Fill(lon, lat, -btheta*1e-5); + hbdown_counter->Fill(lon, lat, -br*1e-5); + hbeast_counter->Fill(lon, lat, bphi*1e-5); + } // This histograms is now filled with the number of entries. // Below we will divide with the time (in seconds) to get // event rate per bin. @@ -555,10 +571,18 @@ TH2F *trigAndOr_rate = (TH2F*) trigAndOr_counter->Clone("trigAndOr_rate"); TH2F *trigAndAnd_rate = (TH2F*) trigAndAnd_counter->Clone("trigAndAnd_rate"); TH2F *nd_rate = (TH2F*) nd_counter->Clone("nd_rate"); - TH2F *hbabs_norm = (TH2F*) hbabs_counter->Clone("hbabs_norm"); - TH2F *hbnorth_norm = (TH2F*) hbnorth_counter->Clone("hbnorth_norm"); - TH2F *hbdown_norm = (TH2F*) hbabs_counter->Clone("hbdown_norm"); - TH2F *hbeast_norm = (TH2F*) hbabs_counter->Clone("hbeast_norm"); + + TH2F *hbabs_norm; + TH2F *hbnorth_norm; + TH2F *hbdown_norm; + TH2F *hbeast_norm; + + if(field) { + hbabs_norm = (TH2F*) hbabs_counter->Clone("hbabs_norm"); + hbnorth_norm = (TH2F*) hbnorth_counter->Clone("hbnorth_norm"); + hbdown_norm = (TH2F*) hbabs_counter->Clone("hbdown_norm"); + hbeast_norm = (TH2F*) hbabs_counter->Clone("hbeast_norm"); + } // Now we divide each histogram _counter with the time histogram // obtBinTime to have an histogram _rate. Note that, when a second @@ -625,26 +649,27 @@ // fill the bins with the values of the magnetic field for each // event, we need to divide with the number of fills done, that is // event_counter. - hbabs_norm->Divide(hbabs_counter, event_counter, 1, 1, ""); - oss.str(""); - oss << basename.Data() << "_orbit_Babs.png"; - printHist(hbabs_norm, mapFile, outDirectory, oss.str().c_str(), "B abs (G)", -width, height, 0, 0); - - hbnorth_norm->Divide(hbnorth_counter, event_counter, 1, 1, ""); - oss.str(""); - oss << basename.Data() << "_orbit_Bnorth.png"; - printHist(hbnorth_norm, mapFile, outDirectory, oss.str().c_str(), "B north (G)", -width, height, 0, 1); - - hbdown_norm->Divide(hbdown_counter, event_counter, 1, 1, ""); - oss.str(""); - oss << basename.Data() << "_orbit_Bdown.png"; - printHist(hbdown_norm, mapFile, outDirectory, oss.str().c_str(), "B down (G)", -width, height, 0, 1); - - hbeast_norm->Divide(hbeast_counter, event_counter, 1, 1, ""); - oss.str(""); - oss << basename.Data() << "_orbit_Beast.png"; - printHist(hbeast_norm, mapFile, outDirectory, oss.str().c_str(), "B east (G)", -width, height, 0, 1); - + if(field) { + hbabs_norm->Divide(hbabs_counter, event_counter, 1, 1, ""); + oss.str(""); + oss << basename.Data() << "_orbit_Babs.png"; + printHist(hbabs_norm, mapFile, outDirectory, oss.str().c_str(), "B abs (G)", -width, height, 0, 0); + + hbnorth_norm->Divide(hbnorth_counter, event_counter, 1, 1, ""); + oss.str(""); + oss << basename.Data() << "_orbit_Bnorth.png"; + printHist(hbnorth_norm, mapFile, outDirectory, oss.str().c_str(), "B north (G)", -width, height, 0, 1); + + hbdown_norm->Divide(hbdown_counter, event_counter, 1, 1, ""); + oss.str(""); + oss << basename.Data() << "_orbit_Bdown.png"; + printHist(hbdown_norm, mapFile, outDirectory, oss.str().c_str(), "B down (G)", -width, height, 0, 1); + + hbeast_norm->Divide(hbeast_counter, event_counter, 1, 1, ""); + oss.str(""); + oss << basename.Data() << "_orbit_Beast.png"; + printHist(hbeast_norm, mapFile, outDirectory, oss.str().c_str(), "B east (G)", -width, height, 0, 1); + } delete obtBinTime; delete event_counter; @@ -668,14 +693,16 @@ delete trigS111A_rate; delete trigS12andS21andS22_rate; - delete hbabs_counter; - delete hbnorth_counter; - delete hbdown_counter; - delete hbeast_counter; - delete hbabs_norm; - delete hbnorth_norm; - delete hbdown_norm; - delete hbeast_norm; + if(field) { + delete hbabs_counter; + delete hbnorth_counter; + delete hbdown_counter; + delete hbeast_counter; + delete hbabs_norm; + delete hbnorth_norm; + delete hbdown_norm; + delete hbeast_norm; + } rootFile->Close(); } @@ -699,7 +726,7 @@ // Scale() and Merge()). // // This function depends on InitStyle(); -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) { InitStyle(); @@ -797,7 +824,7 @@ // querying the database with the RESURS DK-1 id number 29228, // selecting the widest timespan, including the satellite name in the // results. -cTle *getTle(TString tleFile, TDatime offRes) +cTle *getTle(TString tleFile, TTimeStamp offResTS) { Float_t tledatefromfile, tledatefromroot; fstream tlefile(tleFile.Data(), ios::in); @@ -827,7 +854,10 @@ // Sort by date sort(ctles.begin(), ctles.end(), compTLE); - tledatefromroot = (offRes.GetYear()-2000)*1e3 + (offRes.Convert() - (TDatime(offRes.GetYear(), 1, 1, 0, 0, 0)).Convert())/ (24.*3600.); + UInt_t year, month, day; + offResTS.GetDate(kTRUE, 0, &year, &month, &day); + TTimeStamp firstofjan = TTimeStamp(year, 1, 1, 0, 0, 0); + tledatefromroot = (year-2000)*1e3 + (offResTS.GetSec() - firstofjan.GetSec())/(24.*3600.); for(iter = ctles.begin(); iter != ctles.end(); iter++) { cTle *tle = *iter; @@ -871,7 +901,7 @@ // Look for a timesync in the TFile rootFile. Set timesync and // obt_timesync. Returns 1 if timesync is found, 0 otherwise. -int lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) { +UInt_t lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) { *timesync = -1; // will be != -1 if found ULong64_t nevents = 0; @@ -942,8 +972,10 @@ } -// Return a string like YYYY-MM-DD hh:mm:ss, a datetime format. -string getTleDatetime(cTle *tle) +// +// Returns the tle date as a TTimeStamp object. +// +TTimeStamp getTleDatetime(cTle *tle) { int year, mon, day, hh, mm, ss; double dom; // day of month (is double!) @@ -963,9 +995,9 @@ ss = (int) floor( ((((dom - day) * 24 - hh) * 60 - mm) * 60)); // ms = (int) floor( (((((dom - day) * 24 - hh) * 60 - mm) * 60) - ss) * 1000); - date << year << "-" << mon << "-" << day << " " << hh << ":" << mm << ":" << ss; + TTimeStamp t = TTimeStamp(year, mon, day, hh, mm, ss, 0, true); - return date.str(); + return t; } //