--- quicklook/OrbitalRate/src/OrbitalRate.cpp 2006/12/06 15:52:13 1.2 +++ quicklook/OrbitalRate/src/OrbitalRate.cpp 2006/12/07 13:27:30 1.4 @@ -346,19 +346,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 +356,8 @@ 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 + recalc_(y, m*31+d, 0, 0, 0); /********** Magnetic Field **************/ tr = (TTree*)rootFile->Get("Physics"); @@ -469,16 +457,15 @@ 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); + 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<<", "< 1) cout << endl << "******** deltaTime<1 ********" << endl; // Does nothing for the first two events or if acquisition time if more // than 1s. - if(i<2 || (deltaTime > 1)) continue; + if(i<1 || (deltaTime > 1)) continue; // CAS3 and CAS4 are not rates but only counters. So I fill // with the bin with the difference beetween the actual counter @@ -498,12 +485,10 @@ // this values but I need to count how many times I fill // each bin. This is done by the histogram event_counter. // I will normalize later. - hbabs_counter->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); + 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 @@ -574,8 +559,6 @@ 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"); // Now we divide each histogram _counter with the time histogram // obtBinTime to have an histogram _rate. Note that, when a second @@ -591,19 +574,19 @@ trigS111A_rate->Divide(trigS111A_counter, obtBinTime, 1, 1, ""); oss.str(""); oss << basename.Data() << "_orbit_trigS111A.png"; - trigS111A_rate->SetMinimum(10); + trigS111A_rate->SetMinimum(9); printHist(trigS111A_rate, mapFile, outDirectory, oss.str().c_str(), "S111A (Hz)", -width, height, true, 0); antiCAS4_rate->Divide(antiCAS4_counter, obtBinTime, 1, 1, ""); oss.str(""); oss << basename.Data() << "_orbit_CAS4.png"; - antiCAS4_rate->SetMinimum(100); + antiCAS4_rate->SetMinimum(99); printHist(antiCAS4_rate, mapFile, outDirectory, oss.str().c_str(), "CAS4 (Hz)", -width, height, true, 0); antiCAS3_rate->Divide(antiCAS3_counter, obtBinTime, 1, 1, ""); oss.str(""); oss << basename.Data() << "_orbit_CAS3.png"; - antiCAS3_rate->SetMinimum(100); + antiCAS3_rate->SetMinimum(99); printHist(antiCAS3_rate, mapFile, outDirectory, oss.str().c_str(), "CAS3 (Hz)", -width, height, true, 0); event_rate->Divide(event_counter, obtBinTime, 1, 1, ""); @@ -614,13 +597,13 @@ trigS11andS12_rate->Divide(trigS11andS12_counter, obtBinTime, 1, 1, ""); oss.str(""); oss << basename.Data() << "_orbit_trigS11andS12.png"; - antiCAS3_rate->SetMinimum(100); + trigS11andS12_rate->SetMinimum(99); printHist(trigS11andS12_rate, mapFile, outDirectory, oss.str().c_str(), "(S11*S12) (Hz)", -width, height, 1, 0); trigS12andS21andS22_rate->Divide(trigS12andS21andS22_counter, obtBinTime, 1, 1, ""); oss.str(""); oss << basename.Data() << "_orbit_trigS12andS21andS22.png"; - antiCAS3_rate->SetMinimum(10); + trigS12andS21andS22_rate->SetMinimum(9); printHist(trigS12andS21andS22_rate, mapFile, outDirectory, oss.str().c_str(), "(S12*S12*S21) (Hz)", -width, height, true, 0); trigAndOr_rate->Divide(trigAndOr_counter, obtBinTime, 1, 1, ""); @@ -662,16 +645,6 @@ 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); - delete obtBinTime; delete event_counter; @@ -699,14 +672,10 @@ 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; rootFile->Close(); } @@ -737,11 +706,9 @@ // Create a canvas and draw the TH2F with a nice colormap for z // values, using log scale for z values, if requested, and setting // some title. - TCanvas *canvas = new TCanvas("h", "passed histogram", width*2, height*2); + TCanvas *canvas = new TCanvas("h", "h histogram", width*2, height*2); - if(use_log) { - canvas->SetLogz(); - } + if(use_log) canvas->SetLogz(); h->SetTitle(title); h->SetXTitle("Longitude (deg)");