--- quicklook/OrbitalRate/src/OrbitalRate.cpp 2006/12/06 16:25:52 1.3 +++ quicklook/OrbitalRate/src/OrbitalRate.cpp 2007/04/02 19:44:38 1.7 @@ -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; @@ -346,19 +363,7 @@ /********** Magnetic Field **************/ // Check that all this is correct! - float dimo = 0.0; // dipole moment (computed from dat files) - 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_(); + float br, btheta, bphi; // 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 @@ -368,8 +373,9 @@ Int_t d = tledate.GetDay(); float year = (float) y + (m*31+d)/365; - // Compute the magnetic dipole moment. - feldcof_(&year, &dimo); + // Initialize common data for geopack + if(field) + recalc_(y, m*31+d, 0, 0, 0); /********** Magnetic Field **************/ tr = (TTree*)rootFile->Get("Physics"); @@ -416,7 +422,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 @@ -469,9 +475,9 @@ alt = coo.m_Alt; /********** Magnetic Field **************/ - feldg_(&lat, &lon, &alt, &bnorth, &beast, &bdown, &babs); - shellg_(&lat, &lon, &alt, &dimo, &xl, &icode, &bab1); - findb0_(&stps, &bdel, &value, &bequ, &rr0); + 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, babs); - hbnorth_counter->Fill(lon, lat, bnorth); - hbdown_counter->Fill(lon, lat, bdown); - hbeast_counter->Fill(lon, lat, beast); - hb0_counter->Fill(lon, lat, bequ); - hl_counter->Fill(lon, lat, xl); - + 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. @@ -570,12 +575,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 *hb0_norm = (TH2F*) hb0_counter->Clone("hb0_norm"); - TH2F *hl_norm = (TH2F*) hl_counter->Clone("hl_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 @@ -642,36 +653,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); - - 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); - + 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; @@ -695,18 +697,16 @@ delete trigS111A_rate; delete trigS12andS21andS22_rate; - delete hbabs_counter; - delete hbnorth_counter; - delete hbdown_counter; - delete hbeast_counter; - delete hb0_counter; - delete hl_counter; - delete hbabs_norm; - delete hbnorth_norm; - delete hbdown_norm; - delete hbeast_norm; - delete hb0_norm; - delete hl_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(); } @@ -902,7 +902,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;