OrbitalRate/src/OrbitalRate.cpp

/**
 * OrbitalRate
 * author  Nagni
 * version 1.0 - 27 April 2006
 *
 * version 2.0
 * author De Simone 
 * - most of the code rewritten
 * - added graphs, magnetic field, new overflow resolution (AC), tle
 * stuff.
 *
 */
#include <physics/anticounter/AnticounterEvent.h>
#include <physics/trigger/TriggerEvent.h>
#include <physics/neutronDetector/NeutronEvent.h>
#include "physics/neutronDetector/NeutronRecord.h"
#include <mcmd/McmdEvent.h>
#include <mcmd/McmdRecord.h>
#include <EventHeader.h>
#include <PscuHeader.h>
#include <TTree.h>
#include "sgp4.h"
#include "TH2F.h"
#include "TFrame.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TASImage.h"
#include <TDatime.h>
#include <TFile.h>

#include <TTimeStamp.h>
#include "TString.h"
#include "TObjString.h"
#include "TStyle.h"
#include "TPaletteAxis.h"
#include "TROOT.h"
#include <sys/stat.h>
#include <fstream>
#include <iostream>

#include <OrbitalRate.h>

using namespace std;

int main(int argc, char* argv[]){
  TString *rootFile = NULL;
  TString outDir  = "./";
  TString mapFile = "";
  TString tleFile = "";
  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");
    printf("Try '--help' for more information. \n");
    exit(1);
  }  

  if (!strcmp(argv[1], "--help")){
    printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");
    printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");
    printf( "\t --help                 Print this help and exit \n");
    printf( "\t -tle[File path]         Path where to find the tle infos \n");
    printf( "\t\tUse the script retrieve_TLE.sh to create the file.\n ");
    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);
  }

  // Ok, here we should have at least one root file.  We check that
  // the filename contains ".root".
  if(strstr(argv[1], ".root"))
    rootFile = new TString(argv[1]);
  else {
    cerr << "OrbitalRate: no root file." << endl << "See --help" << endl;
    exit(EXIT_FAILURE);
  }

  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");
        printf( "Try '--help' for more information. \n");
        exit(1);
      } else {
        outDir = argv[i];
        continue;
      }
    }

    if (!strcmp(argv[i], "-tle")){
      if (++i >= argc){
        printf( "-tle needs arguments. \n");
        printf( "Try '--help' for more information. \n");
        exit(1);
      } else {
        tleFile = argv[i];
        continue;
      }
    }

    if (!strcmp(argv[i], "-offTime")){
      if (++i >= argc){
        printf( "-offTime needs arguments. \n");
        printf( "Try '--help' for more information. \n");
        exit(1);
      } 
      else{
        offTime = atol(argv[i]);
        continue;
      }
    }

    if (!strcmp(argv[i], "-offDate")){
      if (++i >= argc){
        printf( "-offDate needs arguments. \n");
        printf( "Try '--help' for more information. \n");
        exit(1);
      } 
      else{
        offDate = atol(argv[i]);
        continue;
      }
    }

    if (!strcmp(argv[i], "-map")){
      if (++i >= argc){
        printf( "-map needs arguments. \n");
        printf( "Try '--help' for more information. \n");
        exit(1);
      } else {
        mapFile = argv[i];
        continue;
      }
    }
  }

  if (mapFile != ""){
    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");
  }
}


void InitStyle() {
  gROOT->SetStyle("Plain");

  TStyle *myStyle[2], *tempo;
  myStyle[0]=new TStyle("StyleWhite", "white");
  myStyle[1]=new TStyle("StyleBlack", "black");

  tempo=gStyle;
  Int_t linecol, bkgndcol, histcol;
 
  for(Int_t style=0; style<2; style++) {

    linecol=kWhite*style+kBlack*(1-style);
    bkgndcol=kBlack*style+kWhite*(1-style);
    histcol=kYellow*style+kBlack*(1-style); // was 95

    myStyle[style]->Copy(*tempo);

    myStyle[style]->SetCanvasBorderMode(0);
    myStyle[style]->SetCanvasBorderSize(1);
    myStyle[style]->SetFrameBorderSize(1);
    myStyle[style]->SetFrameBorderMode(0);
    myStyle[style]->SetPadBorderSize(1);
    myStyle[style]->SetStatBorderSize(1);
    myStyle[style]->SetTitleBorderSize(1);
    myStyle[style]->SetPadBorderMode(0);
    myStyle[style]->SetPalette(1,0);
    myStyle[style]->SetPaperSize(20,27);        
    myStyle[style]->SetFuncColor(kRed);
    myStyle[style]->SetFuncWidth(1);
    myStyle[style]->SetLineScalePS(1); 
    myStyle[style]->SetCanvasColor(bkgndcol);
    myStyle[style]->SetAxisColor(linecol,"XYZ");
    myStyle[style]->SetFrameFillColor(bkgndcol);
    myStyle[style]->SetFrameLineColor(linecol);
    myStyle[style]->SetLabelColor(linecol,"XYZ");
    myStyle[style]->SetPadColor(bkgndcol);
    myStyle[style]->SetStatColor(bkgndcol);
    myStyle[style]->SetStatTextColor(linecol);
    myStyle[style]->SetTitleColor(linecol,"XYZ");
    myStyle[style]->SetTitleFillColor(bkgndcol);
    myStyle[style]->SetTitleTextColor(linecol);
    myStyle[style]->SetLineColor(linecol);
    myStyle[style]->SetMarkerColor(histcol);
    myStyle[style]->SetTextColor(linecol);
    
    myStyle[style]->SetGridColor((style)?13:kBlack);
    myStyle[style]->SetHistFillStyle(1001*(1-style));
    myStyle[style]->SetHistLineColor(histcol);
    myStyle[style]->SetHistFillColor((style)?bkgndcol:kYellow);

    myStyle[style]->SetOptStat(0); // Remove statistic summary
  }

  myStyle[1]->cd();

  gROOT->ForceStyle();

}


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;
  //********************************************************/

  TTree                  *tr         = 0;
  TFile *rootFile;
  FILE *f;

  pamela::McmdEvent      *mcmdev     = 0;
  pamela::McmdRecord     *mcmdrc     = 0;
  pamela::EventHeader    *eh         = 0;
  pamela::PscuHeader     *ph       = 0;
  TArrayC                *mcmddata;
  ULong64_t nevents = 0;
  stringstream  oss;

  Float_t timesync = 0, obt_timesync = 0;
  Long64_t offsetTime = 0;
  Long64_t timeElapsedFromTLE = 0;
  Long64_t deltaTime = 0, oldtimeElapsedFromTLE = 0;
  bool a_second_is_over;

  Float_t lon, lat, alt;

  vector<Double_t> vector_trigAndOr;
  vector<Double_t> vector_trigAndAnd;
  vector<Double_t> vector_trigS11andS12;
  vector<Double_t> vector_trigS12andS21andS22;
  vector<Double_t> vector_trigS111A;

  double mean_trigAndOr;
  double mean_trigAndAnd;
  double mean_trigS11andS12;
  double mean_trigS12andS21andS22;
  double mean_trigS111A;

  // We'll use this size for the generated images.
  TImage *tImage=TImage::Open(mapFile);
  int width=(int)(tImage->GetWidth()*0.80);
  int height=(int)(tImage->GetHeight()*0.80);
  delete tImage;

  // This histogram will store time (in seconds) spent in each bin.
  TH2F *obtBinTime = new TH2F("obtBinTime", "Time of acquisition of background data", 360, -180, 180, 180, -90, 90);

  // Now I create histograms longitude x latitude to hold values.  I
  // use the suffix _counter to say that this values are what I read
  // from Pamela and they are not normalized in any way.

  // This historam will store the number of events occurred in each bin.
  TH2F *event_counter = new TH2F("event_counter", "Event rate", 360, -180, 180, 180, -90, 90);
  TH2F *nd_counter = new TH2F("nd_counter", "Upper background neutrons", 360, -180, 180, 180, -90, 90);
  TH2F *antiCAS4_counter = new TH2F("CAS4_counter", "CAS4 rate", 360, -180, 180, 180, -90, 90);
  TH2F *antiCAS3_counter = new TH2F("CAS3_counter", "CAS3 rate", 360, -180, 180, 180, -90, 90);
  TH2F *trigAndOr_counter = new TH2F("trigAndOr_counter", "Rate of triggering in (S11+S12)*(S21+S22)*(S31+S32) configuration", 360, -180, 180, 180, -90, 90);
  TH2F *trigAndAnd_counter = new TH2F("trigAndAnd_counter", "Rate of triggering in (S11*S12)*(S21*S22)*(S31*S32) configuration", 360, -180, 180, 180, -90, 90);
  TH2F *trigS11andS12_counter = new TH2F("trigS11andS12_counter", "Rate of S1 triggers", 360, -180, 180, 180, -90, 90); //(S11+S12)
  TH2F *trigS12andS21andS22_counter = new TH2F("trigS12andS21andS22_counter", "Rate of S11*S21*S21 triggers", 360, -180, 180, 180, -90, 90); //(S11*S12*S21)
  TH2F *trigS111A_counter = new TH2F("trigS111A_counter", "Rate of S111A counts", 360, -180, 180, 180, -90, 90); //(S111A)

  // 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;
  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;
  basename = ((TObjString*) filename->Tokenize('/')->Last())->GetString(); // we get file.root
  basename = ((TObjString*)basename.Tokenize('.')->First())->GetString(); // we get file

  // Exit if the map file doesn't exist.
  if(! (f = fopen(mapFile.Data(), "r")) ) {
    cerr << "Error: the file " << mapFile.Data() << " does not exists." << endl;
    exit(EXIT_FAILURE);
  }

  // Open the root file.
  rootFile = new TFile(filename->Data());
  if (rootFile->IsZombie()) {
    printf("The file %s does not exist\n", (filename->Data()));
    exit(EXIT_FAILURE);
  }

  // Look for a timesync in the TFile rootFile.  We also get the obt
  // of the timesync mcmd.
  bool err;
  err = lookforTimesync(rootFile, &timesync, &obt_timesync);
  if(!err) {
    cerr << "Warning!!! No timesync info has been found in the file "
         << filename->Data() << endl;
    exit(EXIT_FAILURE);
  }

  // Here I do: resurs offset + timesync
  TDatime offRes = TDatime(offDate, offTime);
  TTimeStamp offResTS = TTimeStamp(offRes.GetYear(), offRes.GetMonth(), offRes.GetDay(), offRes.GetHour(), offRes.GetMinute(), offRes.GetSecond(), 0, kTRUE, timesync);
//   cout<<"timesync="<<timesync<<"  obt_timesync="<<obt_timesync<<endl;
//   cout<<"offResTS.GetSec()="<<offResTS.GetSec()<<endl;
//   cout<<"offResTS.AsString()="<<offResTS.AsString()<<endl<<endl;

  // 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.
  string str1 = "RESURS-DK 1";
  string str2 = "1 29228U 06021A   06170.19643714  .00009962  00000-0  21000-3 0   196";
  string str3 = "2 29228 069.9363 054.7893 0167576 127.4359 017.0674 15.31839265   604";
  cTle *tle1 = new cTle(str1, str2, str3);

  // If we have to use a TLE file, call getTle().
  if (tleFile != "")
    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
  //  cout<<"here0 "<<offResTS.GetSec()<<endl;
  //  offResTS.Set(offResTS.GetSec() - obt_timesync, kTRUE, 0, kFALSE);
  offResTS.Set(offResTS.GetSec() - obt_timesync, kFALSE, 0, kFALSE);
  //  cout<<"here1 "<<offResTS.GetSec()<<endl;

  cOrbit       orbit(*tle1);
  cEci         eci;
  cCoordGeo    coo;

  // 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) offResTS.GetSec() - (Long64_t) tledate.GetSec());
//   cerr<<"offsetTime="<<offsetTime<<endl
//       <<"offResTS.GetSec()="<<offResTS.GetSec()<<endl
//       <<"tledate.GetSec()="<<tledate.GetSec()<<endl<<endl;

  /********** Magnetic Field **************/
  // Check that all this is correct!
  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
  // beacause changes are not relevant at all.
  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
  if(field) 
    recalc_(y, m*31+d, 0, 0, 0);
  /********** Magnetic Field **************/

  tr = (TTree*)rootFile->Get("Physics");
  TBranch *headBr = tr->GetBranch("Header");
  tr->SetBranchAddress("Header", &eh);

  /********** Anticounter **************/
  pamela::anticounter::AnticounterEvent *antiev = 0;
  tr->SetBranchAddress("Anticounter", &antiev);

  Int_t oldCAS4 = 0;
  Int_t diffCAS4 = 0;
  Int_t oldCAS3 = 0;
  Int_t diffCAS3 = 0;
  /********** Anticounter **************/
        
  /********** Trigger **************/
  pamela::trigger::TriggerEvent *trigger   = 0;
  tr->SetBranchAddress("Trigger", &trigger);

  Int_t oldtrigAndOr = 0;
  Int_t oldtrigAndAnd = 0;
  Int_t oldtrigS11andS12 = 0;
  Int_t oldtrigS12andS21andS22 = 0;
  Int_t oldtrigS111A = 0;
  /********** Trigger **************/

  /********** ND **************/
  Int_t    tmpSize=0;
  Int_t    sumTrig=0; 
  Int_t    sumUpperBackground=0; 
  Int_t    sumBottomBackground=0;

  pamela::neutron::NeutronRecord *nr = 0;
  pamela::neutron::NeutronEvent *ne = 0;
  tr->SetBranchAddress("Neutron", &ne);
  /********** ND **************/
        
  nevents = tr->GetEntries();

  for(UInt_t i = 0; i < nevents; i++) //Fill variables from root-ple
    { 
      tr->GetEntry(i);
      ph = eh->GetPscuHeader();

      // obt in ms
      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
      // because the indetermination on the timesync is about 1s.
      timeElapsedFromTLE = offsetTime + obt/1000;
      if(!i) cerr<<"1st event: timeElapsedFromTLE="<<timeElapsedFromTLE<<endl;

      // I also need the abstime in seconds rounded to the lower
      // value.  Every second, we set a_second_is_over to true.  Only
      // in this case histograms with triggers are filled.
      a_second_is_over = (timeElapsedFromTLE > oldtimeElapsedFromTLE) ? 1 : 0;
      oldtimeElapsedFromTLE = timeElapsedFromTLE;

      // I need the acquisition time between two triggers to fill the
      // obtBinTime (histo of time spent in the bin).  The time is in
      // second.
      deltaTime = timeElapsedFromTLE - oldtimeElapsedFromTLE;
      oldtimeElapsedFromTLE = timeElapsedFromTLE;

      // Finally, we get coordinates from absolute time the orbit
      // object initialised with the TLE data.  cOrbit::getPosition()
      // requires the elapased time from the tle in minutes.
      // Coordinates are stored in the structure eci.
      orbit.getPosition(((double) timeElapsedFromTLE)/60., &eci);
      coo = eci.toGeo();

      /********** ND **************/
      // Summing over all stored pamela::neutron::NeutronRecords in
      // this event *ne.
      for(Int_t j = 0; j < ne->Records->GetEntries(); j++) {
        nr = (pamela::neutron::NeutronRecord*)ne->Records->At(j);
        sumTrig += (int)nr->trigPhysics;
        sumUpperBackground += (int)nr->upperBack;
        sumBottomBackground += (int)nr->bottomBack;
      }
      /********** ND **************/

      /********** Anticounter **************/
      // Get the difference between the actual counter and the
      // previous counter for anticoincidence, dealing with the
      // overflow with solve_ac_overflow().
      diffCAS4 = solve_ac_overflow(oldCAS4, antiev->counters[0][6]);
      diffCAS3 = solve_ac_overflow(oldCAS3, antiev->counters[0][10]);
      /********** Anticounter **************/

      // Build coordinates in the right range.  We want to convert,
      // just for aesthetic, longitude from (0, 2*pi) to (-pi, pi).
      // We also want to convert from radians to degrees.
      lon = (coo.m_Lon > PI) ? rad2deg(coo.m_Lon - 2*PI) : rad2deg(coo.m_Lon);
      lat = rad2deg(coo.m_Lat);
      alt = coo.m_Alt;

      /********** Magnetic Field **************/
      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<<", "<<M_PI/2.-coo.m_Lat<<", "<<coo.m_Lon<<")"<<endl;
      /********** Magnetic Field **************/

      // serve fare il controllo deltatime < 1?
      if (deltaTime > 1) cout << endl << "******** deltaTime<1 ********" << endl;
      // Does nothing for the first two events or if acquisition time if more
      // than 1s.
      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
      // and the previous one and then divide with the time (see
      // below) to have rates.
      if(diffCAS3>1e3) // additional cut to avoid the peaks after dead time
        diffCAS3 = (Int_t) antiCAS3_counter->GetBinContent((Int_t)antiCAS3_counter->GetEntries()-1);
      antiCAS3_counter->Fill(lon , lat, diffCAS3);

      if(diffCAS4>1e3) // additional cut to avoid the peaks after dead time
        diffCAS4 = (Int_t) antiCAS4_counter->GetBinContent((Int_t) antiCAS4_counter->GetEntries()-1); 
      antiCAS4_counter->Fill(lon, lat, diffCAS4);

      // Magnetic field values should be handled a bit carefully.
      // For every event I get a position and the related magnetic
      // field values.  I can fill the histograms lon x lat with
      // 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.
      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.
      event_counter->Fill(lon, lat);

      // counters about triggers are already rates (Hz).  Only
      // every second we fill fill with the mean over all values.
      if(a_second_is_over) {
        // This histograms will hold the time, in seconds, spent
        // in the bin.
        obtBinTime->Fill(lon, lat, 1);

        // get the means
        mean_trigAndOr = getMean(vector_trigAndOr);
        mean_trigAndAnd = getMean(vector_trigAndAnd);
        mean_trigS11andS12 = getMean(vector_trigS11andS12);
        mean_trigS12andS21andS22 = getMean(vector_trigS12andS21andS22);
        mean_trigS111A = getMean(vector_trigS111A);

        // clear data about the last second
        vector_trigAndOr.clear();
        vector_trigAndAnd.clear();
        vector_trigS11andS12.clear();
        vector_trigS12andS21andS22.clear();
        vector_trigS111A.clear();

        // Fill with the mean rate value
        trigAndOr_counter->Fill(lon , lat, mean_trigAndOr);
        trigAndAnd_counter->Fill(lon , lat, mean_trigAndAnd);
        trigS11andS12_counter->Fill(lon , lat, mean_trigS11andS12);
        trigS12andS21andS22_counter->Fill(lon , lat, mean_trigS12andS21andS22);
        trigS111A_counter->Fill(lon, lat, mean_trigS111A);
      }
      else { // Collect values for all the second
        vector_trigAndOr.push_back((1/4.)*trigger->trigrate[0]);
        vector_trigAndAnd.push_back((1/4.)*trigger->trigrate[1]);
        // pmtpl[0] is the rate every 60ms but I want Hz.
        vector_trigS11andS12.push_back((1000./60.)*trigger->pmtpl[0]);
        vector_trigS12andS21andS22.push_back((1/4.)*trigger->trigrate[4]);
        vector_trigS111A.push_back(1.*trigger->pmtcount1[0]);
      }

      // Now we discard ND data if:
      // - NeutronEvent is corrupted.
      if((ne->unpackError != 1))
        nd_counter->Fill(lon, lat, 1.*(sumUpperBackground+sumTrig));

      // Reset counters for ND.
      sumTrig = 0;
      sumUpperBackground = 0;
      sumBottomBackground = 0;
    }

  // We now need to normalize the histograms to print something
  // meaningful.  I create similar histograms with the suffix _rate or
  // _norm.
  TH2F *event_rate = (TH2F*) event_counter->Clone("event_rate");
  TH2F *trigS111A_rate = (TH2F*) trigS111A_counter->Clone("trigS111A_rate");
  TH2F *antiCAS4_rate = (TH2F*) antiCAS4_counter->Clone("antiCAS4_rate");
  TH2F *antiCAS3_rate = (TH2F*) antiCAS3_counter->Clone("antiCAS3_rate");
  TH2F *trigS11andS12_rate = (TH2F*) trigS11andS12_counter->Clone("trigS11andS12_rate");
  TH2F *trigS12andS21andS22_rate = (TH2F*) trigS12andS21andS22_counter->Clone("trigS12andS21andS22_rate");
  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 *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
  // is passed in the above cycle, we fill the histogram obtBinTime
  // with 1 (second) together with all the other histograms.  So
  // dividing here does make sense.
  //
  // Then we call printHist() for each filled TH2F.  These are
  // refered to the root file we're now reading.  We also build up a
  // filename to be passed to the function.  Pay attention that the
  // filename must end with a file format (such as .png or .pdf)
  // recognised by TPad::SaveAs().
  trigS111A_rate->Divide(trigS111A_counter, obtBinTime, 1, 1, "");
  oss.str("");
  oss << basename.Data() << "_orbit_trigS111A.png";
  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(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(99);
  printHist(antiCAS3_rate, mapFile, outDirectory, oss.str().c_str(), "CAS3 (Hz)", -width, height, true, 0);

  event_rate->Divide(event_counter, obtBinTime, 1, 1, "");
  oss.str("");
  oss << basename.Data() << "_orbit_EventRate.png";
  printHist(event_rate, mapFile, outDirectory, oss.str().c_str(), "Event rate (Hz)", -width, height, 0, 0);

  trigS11andS12_rate->Divide(trigS11andS12_counter, obtBinTime, 1, 1, "");
  oss.str("");
  oss << basename.Data() << "_orbit_trigS11andS12.png";
  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";
  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, "");
  oss.str("");
  oss << basename.Data() << "_orbit_trigANDofOR.png";
  printHist(trigAndOr_rate, mapFile, outDirectory, oss.str().c_str(), "(S11+S12)*(S21+S22)*(S31+S32) (Hz)", -width, height, 0, 0);

  trigAndAnd_rate->Divide(trigAndAnd_counter, obtBinTime, 1, 1, "");
  oss.str("");
  oss << basename.Data() << "_orbit_trigANDofAND.png";
  printHist(trigAndAnd_rate, mapFile, outDirectory, oss.str().c_str(), "(S11*S12)*(S21*S22)*(S31*S32) (Hz)", -width, height, 0, 0);

  nd_rate->Divide(nd_counter, obtBinTime, 1, 1, "");
  oss.str("");
  oss << basename.Data() << "_orbit_ND.png";
  printHist(nd_rate, mapFile, outDirectory, oss.str().c_str(), "Neutron rate (Hz)", -width, height, 0, 0);

  // Also normalize histograms about magnetic fields.  Beacause we
  // 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.
  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;

  delete nd_counter;
  delete antiCAS4_counter;
  delete antiCAS3_counter;
  delete trigAndOr_counter;
  delete trigAndAnd_counter;
  delete trigS11andS12_counter;
  delete trigS111A_counter;
  delete trigS12andS21andS22_counter;

  delete event_rate;
  delete nd_rate;
  delete antiCAS4_rate;
  delete antiCAS3_rate;
  delete trigAndOr_rate;
  delete trigAndAnd_rate;
  delete trigS11andS12_rate;
  delete trigS111A_rate;
  delete trigS12andS21andS22_rate;

  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();
}


// Print the istogram *h on the file outputfilename in the direcotry
// outDirectory, using mapFile as background image, sizing the image
// width per height.  Log scale will be used if use_log is true.
//
// If bool_shift is true a further process is performed to solve a
// problem with actual root version (5.12).  This should be used when
// the histrogram is filled also with negative values, because root
// draws zero filled bins (so I have all the pad colorized and this is
// really weird!).  To avoid this problem I shift all the values in a
// positive range and draw again using colz.  Now I will not have zero
// filled bins painted but the scale will be wrong.  This is why I
// need to draw a new axis along the palette.
//
// Pay attention: you cannot use a mapFile different from the provided
// one without adjusting the scaling and position of the image (see
// Scale() and Merge()).
//
// This function depends on InitStyle();
int printHist(TH2F *h, TString mapFile, TString outDirectory, TString outputFilename, const char *title, int width, int height, bool use_log, bool bool_shift)
{
  InitStyle();

  // 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", "h histogram", width*2, height*2);

  if(use_log) canvas->SetLogz();

  h->SetTitle(title);
  h->SetXTitle("Longitude (deg)");
  h->SetYTitle("Latitude (deg)");
  h->SetLabelColor(0, "X");
  h->SetAxisColor(0, "X");
  h->SetLabelColor(0, "Y");
  h->SetAxisColor(0, "Y");
  h->SetLabelColor(0, "Z");
  h->SetAxisColor(0, "Z");

  h->Draw("colz");
  canvas->Update(); // Update!  Otherwise we can't get any palette.

  // If shift in a positive range required (see comment above).
  if(bool_shift) {
    // Remember the minimum and maximum in this graph.
    Float_t min = h->GetMinimum();
    Float_t max = h->GetMaximum();

    // Shift the graph up by 100.  Increase the value if you still get
    // negative filled bins.
    h = shiftHist(h, 100.0);
    h->SetMinimum(min+100.0);
    h->SetMaximum(max+100.0);
 
    // Hide the current axis of the palette
    TPaletteAxis *palette = (TPaletteAxis*) h->GetListOfFunctions()->FindObject("palette");
    if(!palette) cout << "palette is null" << endl;
    TGaxis *ax = (TGaxis*) palette->GetAxis();
    if(!ax) cout << "ax is null" << endl;
    ax->SetLabelOffset(999);
    ax->SetTickSize(0);

    // Create a new axis of the palette using the right min and max and draw it.
    TGaxis *gaxis = new TGaxis(palette->GetX2(), palette->GetY1(), palette->GetX2(), palette->GetY2(), min, max, 510,"+L");
    gaxis->SetLabelColor(0);
    gaxis->Draw();

    // Update again.
    canvas->Update();
  }

  // We merge two images: the image of the earth read from a file on
  // that one of the TPad of canvas (the histogram).  The first one is
  // scaled and adjusted to fit well inside the frame of the second
  // one.  Finally we draw them both.
  //
  // Here there's a trick to avoid blurring during the merging
  // operation.  We get the image from a canvas sized (width*2 x
  // height*2) and draw it on a canvas sized (width x height).

  TCanvas *mergeCanvas = new TCanvas("", "", width, height);
  TImage *img = TImage::Create();
  TImage *terra = TImage::Create();
  img->FromPad(canvas); // get the TCanvas canvas as TImage
  terra->ReadImage(mapFile, TImage::kPng); // get the png file as TImage
  terra->Scale(1304,830);
  img->Merge(terra, "add", 166, 112); // add image terra to image img
  img->Draw("X"); // see what we get, eXpanding img all over mergeCanvas.

  stringstream oss;
  oss << outDirectory.Data() << "/" << outputFilename.Data();
 
  mergeCanvas->SaveAs(oss.str().c_str());
  mergeCanvas->Close();
  canvas->Close();

  return EXIT_SUCCESS;
}

void saveHist(TH1 *h, TString savetorootfile)
{
  TFile *file = new TFile(savetorootfile.Data(), "update");

  h->Write();
  file->Close();
}


// Get the TLE from tleFile.  The right TLE is that one with the
// closest previous date to offRes, that is the date at the time of
// the first timesync found in the root file.
//
// Warning: you must use a tle file obtained by space-track.org
// 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, TTimeStamp offResTS)
{
  Float_t tledatefromfile, tledatefromroot;
  fstream tlefile(tleFile.Data(), ios::in);
  vector<cTle*> ctles;
  vector<cTle*>::iterator iter;


  // Build a vector of *cTle  
  while(1) {
    cTle *tlef;
    string str1, str2, str3;

    getline(tlefile, str1);
    if(tlefile.eof()) break; 

    getline(tlefile, str2);
    if(tlefile.eof()) break;

    getline(tlefile, str3);
    if(tlefile.eof()) break;

    // We now have three good lines for a cTle.
    tlef = new cTle(str1, str2, str3);
    ctles.push_back(tlef);
  }

  // Sort by date
  sort(ctles.begin(), ctles.end(), compTLE);

  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;

    tledatefromfile = getTleJulian(tle);

    if(tledatefromroot > tledatefromfile) {
      tlefile.close();
      cTle *thisTle = tle;
      ctles.clear();

      return thisTle;
    }
  }

  // File ended withoud founding a TLE with a date after offRes.  We'll use the last aviable date.
  cerr << "Warning: using last available TLE in " << tleFile.Data() << ".  Consider updating your tle file." << endl;

  tlefile.close();
  cTle *thisTle = ctles[ctles.size()-1];
  ctles.clear();

  return thisTle; 
}


// Return whether the first TLE is older than the second
bool compTLE (cTle *tle1, cTle *tle2)
{
  return getTleJulian(tle1) > getTleJulian(tle2);
}


// Return the date of the tle using the format (year-2000)*1e3 +
// julian day.  e.g. 6364.5 is the 31th Dec 2006 12:00:00.
// It does *not* return a cJulian date.
float getTleJulian(cTle *tle) {
  return tle->getField(cTle::FLD_EPOCHYEAR)*1e3 + tle->getField(cTle::FLD_EPOCHDAY);
}


// Look for a timesync in the TFile rootFile.  Set timesync and
// obt_timesync.  Returns 1 if timesync is found, 0 otherwise.
UInt_t lookforTimesync(TFile *rootFile, Float_t *timesync, Float_t *obt_timesync) {
  *timesync = -1;  // will be != -1 if found

  ULong64_t             nevents    = 0;
  pamela::McmdRecord     *mcmdrc     = 0;
  pamela::McmdEvent      *mcmdev     = 0;
  TArrayC                *mcmddata;
  TTree *tr = (TTree*) rootFile->Get("Mcmd");

  tr->SetBranchAddress("Mcmd", &mcmdev);

  nevents = tr->GetEntries();

  // Looking for a timesync.  We stop at the first one found.
  long int recEntries;

  for(UInt_t i = 0; i < nevents; i++) {
    tr->GetEntry(i);
    recEntries = mcmdev->Records->GetEntries();

    for(UInt_t j = 0; j < recEntries; j++) {
      mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);

      if ((mcmdrc != 0) && (mcmdrc->ID1 == 0xE0)) //Is it a TimeSync?
        {
          mcmddata = mcmdrc->McmdData;
          *timesync = (((unsigned int)mcmddata->At(0)<<24)&0xFF000000)
            + (((unsigned int)mcmddata->At(1)<<16)&0x00FF0000) 
            + (((unsigned int)mcmddata->At(2)<<8)&0x0000FF00) 
            + (((unsigned int)mcmddata->At(3))&0x000000FF);
          *obt_timesync = (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);

          goto out;  // a timesync is found
        }
    }
  }

 out:
        
  if (*timesync == -1)
    return 0;
  else
    return 1;
}


// Returns the mean value of the elements stored in the vector v.
double getMean(vector<Double_t> v)
{
  double mean = 0;

  for(int i=0; i < v.size(); i++)
    mean += v.at(i);

  return mean/v.size();
}


// Shift all non zero bins by shift.
TH2F* shiftHist(TH2F* h, Float_t shift)
{
  // Global bin number.
  Int_t nBins = h->GetBin(h->GetNbinsX(), h->GetNbinsY());

  for(int i = 0; i < nBins; i++)
    if(h->GetBinContent(i)) h->AddBinContent(i, shift);

  return h;
}


//
// 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!)
  stringstream date; // date in datetime format

  // create a cJulian from the date in tle
  cJulian jdate = cJulian( 2000 + (int) tle->getField(cTle::FLD_EPOCHYEAR), tle->getField(cTle::FLD_EPOCHDAY)); 

  // get year, month, day of month
  jdate.getComponent(&year, &mon, &dom);

  // build a datetime YYYY-MM-DD hh:mm:ss
  date.str("");
  day = (int) floor(dom);
  hh = (int) floor( (dom - day) * 24);
  mm = (int) floor( ((dom - day) * 24 - hh) * 60);
  ss = (int) floor( ((((dom - day) * 24 - hh) * 60 - mm) * 60));
  //  ms = (int) floor( (((((dom - day) * 24 - hh) * 60 - mm) * 60) - ss) * 1000);

  TTimeStamp t = TTimeStamp(year, mon, day, hh, mm, ss, 0, true);

  return t;
}

//
// Solve the overflow for anticoincidence because this counter is
// stored in 2 bytes so counts from 0 to 65535.
//
// counter is the actual value.
// oldValue is meant to be the previous value of counter.
//
// Example:
// for(...) {
//   ...
//   corrected_diff = solve_ac_overflow(oldValueOfTheCounter, actualValue);
//   ...
// }
//
//
// Returns the corrected difference between counter and oldValue and
// set oldValue to the value of counter.
// Attention: oldValue is a reference.
Int_t solve_ac_overflow(Int_t& oldValue, Int_t counter)
{
  Int_t truediff = 0;

  if (counter < oldValue)  // overflow!
    truediff = 0xFFFF - oldValue + counter;
  else
    truediff = counter - oldValue;

  oldValue = counter;

  return truediff;
}