--- quicklook/OrbitalRate/src/OrbitalRate.cpp 2006/12/23 21:34:58 1.5 +++ 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,7 +215,7 @@ } -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; @@ -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; @@ -357,7 +374,8 @@ 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 +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 @@ -457,7 +475,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 +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 *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 +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); - + 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 +697,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(); } @@ -871,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;