//------------------------------------------------------------------------//
//                    getpedTripT.cxx by R. Sacco                         //
//                    Improvements  by I. Danko                           //
//------------------------------------------------------------------------//
#include <iostream>
#include <fstream>
#include <string>
//#include <boost/regex.hpp>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <sstream>

#include "TGraph.h"
#include "TMultiGraph.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TFile.h"
//#include <ITFBChannelId.hxx> 

int main(int argc, char* argv[]){
gROOT->SetStyle("Plain");

 gStyle->SetHistLineWidth(1);
 gStyle->SetPalette(1);
 gStyle->SetPaperSize(20,26);
 
 // Turn off some borders
 gStyle->SetCanvasBorderMode(0);
 gStyle->SetFrameBorderMode(0);
 gStyle->SetPadBorderMode(0);
 gStyle->SetDrawBorder(0);
 gStyle->SetCanvasBorderSize(0);
 gStyle->SetFrameBorderSize(0);
 gStyle->SetPadBorderSize(1);
 gStyle->SetTitleBorderSize(0);
 
 // Say it in black and white!
 gStyle->SetAxisColor(1, "xyz");
 gStyle->SetCanvasColor(0);
 gStyle->SetFrameFillColor(0);
 gStyle->SetFrameLineColor(1);
 gStyle->SetHistFillColor(0);
 gStyle->SetHistLineColor(1);
 //gStyle->SetPadColor(1);
 gStyle->SetPadColor(kWhite);
 gStyle->SetStatColor(0);
 gStyle->SetStatTextColor(1);
 gStyle->SetTitleColor(1);
 gStyle->SetTitleTextColor(1);
 gStyle->SetLabelColor(1,"xyz");
 gStyle->SetLabelOffset(0.02);
 // Show functions in red...
 gStyle->SetFuncColor(2);

 std::cout << "Number of command line arguments: " << argc - 1 << std::endl;
 //for( int i = 1; i < argc; i++ ) 
 //  std::cout << " argc[" << i << "] = " << argv[i] << std::endl;

 if( argc - 1 < 2 ){
   std::cout << "Not enough arguments! " << std::endl;
   std::cout << "Usage: CalibOffline_getpedTripT.exe <detector> <file list> [<ref file>]" << std::endl;
   std::cout << "  <detector>:  ECAL/P0D/SMRD" << std::endl;
   std::cout << "  <file list>: path to the list of constants files" << std::endl;
   std::cout << "  <ref file>:  path to the reference constants file" << std::endl;
   std::cout << "               (optional, if not specified the first file from the list will be used)" << std::endl;
   return -1;
 }
 
 // 1st argument: detector name
 std::string det = (std::string)argv[1];
 std::cout << "Detector: " << det << std::endl;

 // 2nd argument: constants file list 
 std::string ListName = (std::string)argv[2];
 std::cout << "Input list: " << ListName << std::endl;

 // Name of the reference constants file:
 // 3rd argument if specified at the command line,
 // otherwise it is the the first file from the file list (see later)
 char RefName[1000];
 if( argc > 3 ) strcpy(RefName, argv[3]);


 unsigned long firstCh = 0;
 unsigned short int nRMM = 0;
 unsigned int NunusedCh = 0;
 int ntfb[12];
 unsigned int MaxPed = 20000;
 unsigned int MinPed = 12000;
 float MaxDrift = 250.;
 float Ylow = -250.;
 float Yhigh= 250.;

 if( det == "P0D" ){
    //firstCh = ITFBChannelId(1,0,0,0,0).AsUInt();
    firstCh = 2181041152u;
    nRMM = 6;
    NunusedCh = 736;
    ntfb[0]  = 45;
    ntfb[1]  = 29;
    ntfb[2]  = 45;  
    ntfb[3]  = 29;  
    ntfb[4]  = 29;  
    ntfb[5]  = 29;  
    ntfb[6]  = 0;  
    ntfb[7]  = 0;  
    ntfb[8]  = 0;  
    ntfb[9]  = 0;  
    ntfb[10] = 0;  
    ntfb[11] = 0;  
  } else if( det == "ECAL" ){
    //firstCh = ITFBChannelId(4,0,0,0,0).AsUInt();
    firstCh = 2281704448u;
    nRMM = 12;
    NunusedCh = 1107;
    ntfb[0]  = 28;
    ntfb[1]  = 28;
    ntfb[2]  = 15;  
    ntfb[3]  = 44;  
    ntfb[4]  = 44;  
    ntfb[5]  = 26;  
    ntfb[6]  = 26;  
    ntfb[7]  = 15;  
    ntfb[8]  = 26;  
    ntfb[9]  = 26;  
    ntfb[10] = 44;  
    ntfb[11] = 44;  
  } else if( det == "SMRD" ){
    //firstCh = ITFBChannelId(7,0,0,0,0).AsUInt();
    firstCh = 2382367744u;
    nRMM = 4;
    NunusedCh = 4654;  
    MaxDrift = 500.;
    Ylow  = -500.;
    Yhigh = 500.;
    ntfb[0]  = 32;
    ntfb[1]  = 32;
    ntfb[2]  = 32;  
    ntfb[3]  = 32;  
    ntfb[4]  = 0;  
    ntfb[5]  = 0;  
    ntfb[6]  = 0;  
    ntfb[7]  = 0;  
    ntfb[8]  = 0;  
    ntfb[9]  = 0;  
    ntfb[10] = 0;  
    ntfb[11] = 0;  
  } else {
    std::cout << "Unknown detector: " << det << std::endl;
    return -1;
  }
  
 //Run over all files in list to get time steps
 std::ifstream ListFile( ListName.c_str() );
 int counter = 0;
 double minTime=2000000000;
 double maxTime=0;
 time_t times[1000] = {0};

 while(true){
   //	counter++;
   //Get a filename from the list; break if we are at EOF
   char buf[1000];
   ListFile>>buf;
   if(!ListFile.good()) break;

   // If the reference file is not specified at the command line
   // use the first file in the list as the reference
   if( argc == 3 && counter == 0 ) strcpy(RefName, buf);

   std::ifstream dataFile(buf);
      
   //Get to validity range string
   dataFile.getline(buf,1000,' ');
   dataFile.getline(buf,1000,' ');
   dataFile.getline(buf,1000,'\'');
   dataFile.getline(buf,1000,'\'');
  
   //Convert to a unixtime (yuck)
   struct tm timeinfo;
   
   char tmBuf[3]="**";
   
   tmBuf[0]=buf[2];tmBuf[1]=buf[3];
   timeinfo.tm_year=atoi(tmBuf)+100;
   tmBuf[0]=buf[5];tmBuf[1]=buf[6];
   timeinfo.tm_mon=atoi(tmBuf)-1;
   tmBuf[0]=buf[8];tmBuf[1]=buf[9];
   timeinfo.tm_mday=atoi(tmBuf);
   tmBuf[0]=buf[11];tmBuf[1]=buf[12];
   timeinfo.tm_hour=atoi(tmBuf);
   tmBuf[0]=buf[14];tmBuf[1]=buf[15];
   timeinfo.tm_min=atoi(tmBuf);
   tmBuf[0]=buf[17];tmBuf[1]=buf[18];
   timeinfo.tm_sec=atoi(tmBuf);
   
   time_t unixTime=mktime(&timeinfo);

   if(unixTime>maxTime)maxTime=unixTime;
   if(unixTime<minTime)minTime=unixTime;
   
   //   std::cout << "UTime 1: "  <<  unixTime << std::endl;
   times[counter]=unixTime;
   //   std::cout << "Time at point " << counter << "  is " << times[counter] << std::endl;
   counter++;
 }
  
 //*****Read disconnected channel list from file
 char * soffroot;
 soffroot = getenv("CALIBOFFLINEROOT");
 std::ifstream unusedFile( Form("%s/share/data/%s_UnusedChannels.txt",soffroot,det.c_str()) );
 unsigned long int unusedCh[NunusedCh];
 if( unusedFile.is_open() ){
   for( unsigned int ch = 0; ch < NunusedCh; ch++){
     unusedFile >> unusedCh[ch];
     // std::cout << ch << " " << unusedCh[ch] << std::endl;
   }
   unusedFile.close();
 } else {
    std::cout << "Cannot open file with unused channel list! " << std::endl;
    return -1;
 }

 //Set histogram x axes to comfortably fit all data
 //with 1 bin per hour
 double histTimeMax = maxTime+(double(maxTime-minTime))*0.05;
 double histTimeMin = minTime-(double(maxTime-minTime))*0.05;
 int histBins = (int)(histTimeMax-histTimeMin)/3600;

 // Vectors for TGraphs
 int subtimes[counter];
 int hioverchsubrun[counter];
 int hiunderchsubrun[counter];
 int looverchsubrun[counter];
 int lounderchsubrun[counter];

   
 //Load reference pedestals for each channel and store them into an array
 float ref_hiped[nRMM][48][4][16];
 float ref_loped[nRMM][48][4][16];

 std::cout << "Reference file: " << RefName << std::endl;
 std::ifstream refFile( RefName );
   
 int skipfirstline = 0;
 
 while(refFile.good()){
   char buf[1000];
   refFile.getline(buf,1000);
   
   if (skipfirstline) {
     std::istringstream line(buf);
     char *token;
     long int chan;
     long int channel;
     int rmm=-99, tfb=-99, tript=-99, chid=-99;
     int count=0;
     
     token = strtok( buf," ");
     while( token != NULL )
       {
   /* While there are tokens in "string" */
   if(count==0) {
     chan=strtoul(token,NULL,0);
     channel = chan - firstCh;
     rmm = int(channel/524288);
     tfb = int((channel-rmm*524288)/4096);
     tript = int((channel-rmm*524288-tfb*4096)/32);
     chid = (channel-rmm*524288-tfb*4096)-tript*32;
   } else if (count <= 16) {
     ref_hiped[rmm][tfb][tript][chid + count - 1] = strtof(token,NULL);
   } else if (count <= 32) {
     ref_loped[rmm][tfb][tript][chid + count - 17]=strtof(token,NULL);
   } else {
     std::cout << "Unrecognized case..." << std::endl;
   }
   
   /* Get next token: */
   token = strtok( NULL, " " );
   count++;
       }
     
   } // end if(skipfirstline)
   
   skipfirstline = 1;
   
 }
  
 //Output file for all plots
 // TFile fOut("DQPEDPlots.root","RECREATE");
 TFile fOut(Form("DQPedPlots%s.root",det.c_str()),"RECREATE");

 //One bin per hour and per 0.01 ADC
 int Ybins = (int)(Yhigh - Ylow);

 char LhiName[20]="LoPedDriftAll";
 TH2F* dptLoPedDriftPlot=new TH2F(LhiName,"Low Ped Drift All",histBins,histTimeMin,histTimeMax,Ybins,Ylow,Yhigh);
 dptLoPedDriftPlot->GetXaxis()->SetTimeDisplay(1);
 dptLoPedDriftPlot->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
 dptLoPedDriftPlot->GetXaxis()->SetTimeOffset(0,"gmt");

 char HhiName[20]="HiPedDriftAll";
 TH2F* dptHiPedDriftPlot=new TH2F(HhiName,"High Ped Drift All",histBins,histTimeMin,histTimeMax,Ybins,Ylow,Yhigh);
 dptHiPedDriftPlot->GetXaxis()->SetTimeDisplay(1);
 //dptHiPedDriftPlot->GetXaxis()->SetLabelOffset(0.02);
 dptHiPedDriftPlot->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
 dptHiPedDriftPlot->GetXaxis()->SetTimeOffset(0,"gmt");

 TH2F* LoPedAll = new TH2F("LoPedAll","Low Ped All",histBins,histTimeMin,histTimeMax,200,5000.,25000.);
 LoPedAll->GetXaxis()->SetTimeDisplay(1);
 LoPedAll->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
 LoPedAll->GetXaxis()->SetTimeOffset(0,"gmt");

 TH2F* HiPedAll = new TH2F("HiPedAll","High Ped All",histBins,histTimeMin,histTimeMax,200,5000.,25000.);
 HiPedAll->GetXaxis()->SetTimeDisplay(1);
 HiPedAll->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
 HiPedAll->GetXaxis()->SetTimeOffset(0,"gmt");

 TH2F* dptLoPedDriftPlotRMM[nRMM];
 TH2F* dptHiPedDriftPlotRMM[nRMM];
 TH2F* dptLoPedPlotRMM[nRMM];
 TH2F* dptHiPedPlotRMM[nRMM];
 for( unsigned short int r = 0; r < nRMM; r++){
   std::string lodriftName = Form("LoPedDriftRMM%d",r);
   dptLoPedDriftPlotRMM[r] = new TH2F(lodriftName.c_str(),lodriftName.c_str(),histBins,histTimeMin,histTimeMax,Ybins,Ylow,Yhigh);
   dptLoPedDriftPlotRMM[r]->GetXaxis()->SetTimeDisplay(1);
   dptLoPedDriftPlotRMM[r]->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
   dptLoPedDriftPlotRMM[r]->GetXaxis()->SetTimeOffset(0,"gmt");

   std::string hidriftName = Form("HiPedDriftRMM%d",r);
   dptHiPedDriftPlotRMM[r] = new TH2F(hidriftName.c_str(),hidriftName.c_str(),histBins,histTimeMin,histTimeMax,Ybins,Ylow,Yhigh);
   dptHiPedDriftPlotRMM[r]->GetXaxis()->SetTimeDisplay(1);
   dptHiPedDriftPlotRMM[r]->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
   dptHiPedDriftPlotRMM[r]->GetXaxis()->SetTimeOffset(0,"gmt");

   std::string loname = Form("LoGainPedsRMM%d",r);
   dptLoPedPlotRMM[r] = new TH2F(loname.c_str(),loname.c_str(),histBins,histTimeMin,histTimeMax,300,0.,30000);
   dptLoPedPlotRMM[r]->GetXaxis()->SetTimeDisplay(1);
   dptLoPedPlotRMM[r]->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
   dptLoPedPlotRMM[r]->GetXaxis()->SetTimeOffset(0,"gmt");

   std::string hiname = Form("HiGainPedsRMM%d",r);
   dptHiPedPlotRMM[r] = new TH2F(hiname.c_str(),hiname.c_str(),histBins,histTimeMin,histTimeMax,300,0.,30000);
   dptHiPedPlotRMM[r]->GetXaxis()->SetTimeDisplay(1);
   dptHiPedPlotRMM[r]->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
   dptHiPedPlotRMM[r]->GetXaxis()->SetTimeOffset(0,"gmt");
 }
 
  //std::string newListName = Form("%c_pedlist.dat",det[0]);
 std::ifstream newListFile(ListName.c_str());
 
 int pedcounter=0;
 
 float hiped[nRMM][48][4][16];
 float loped[nRMM][48][4][16];
 
 while(true){
   //	counter++;
   //Get a filename from the list; break if we are at EOF
   char buf[1000];
   newListFile>>buf;
   if(!newListFile.good()) break;
   std::ifstream dataFile(buf);
   std::cout << " File: "  << buf << std::endl;
      
   //	dataFile.getline(buf,1000);
      
   int skipfirstline = 0;

   int hioverflow = 0;
   int hiunderflow = 0;
   int looverflow = 0;
   int lounderflow = 0;

   for (int rmm=0; rmm<nRMM; rmm++){
     for(int tfb = 0; tfb<48; tfb++){
       for(int tript = 0; tript<=3; tript++){
   for(int chid = 0; chid<16; chid++){
     hiped[rmm][tfb][tript][chid] = 0.0;
     loped[rmm][tfb][tript][chid] = 0.0;
   }
       }
     }
   }


   while(dataFile.good()){  

     dataFile.getline(buf,1000);
     
     if (skipfirstline) {	  
       std::istringstream line(buf); 
       //       std::cout << buf << std::endl;
       char *token;
       long int chan;
       long int channel;
       int rmm=-99, tfb=-99, tript=-99, chid=-99;
       //	  float loped,hiped;
       int count=0;
       token = strtok( buf," ");
       while( token != NULL )
   {
     /* While there are tokens in "string" */
     if(count==0) {
       chan=strtoul(token,NULL,0);
       //		  std::cout << "chan = " << chan << std::endl; 
       channel = chan - firstCh;
       rmm = int(channel/524288);
       tfb = int((channel-rmm*524288)/4096);
       tript = int((channel-rmm*524288-tfb*4096)/32);
       chid = (channel-rmm*524288-tfb*4096)-tript*32;
     } else if (count <= 16) {
       hiped[rmm][tfb][tript][chid + count - 1] = strtof(token,NULL);
     } else if (count <= 32) {
       loped[rmm][tfb][tript][chid + count - 17]=strtof(token,NULL);
     } else {
       std::cout << "Unrecognized case..." << std::endl;
     }
     /* Get next token: */
     token = strtok( NULL, " " );
     count++;
   }
     } // end of if(skipfirstline)
     
     skipfirstline = 1;
     
   } //end of while(dataFile.good())
   
   unsigned long int chan = 0;

   for (int rmm=0; rmm<nRMM; rmm++){
     for(int tfb = 0; tfb<ntfb[rmm]; tfb++){
       //***** Skip empty RMM ports in P0D:
       if( det == "P0D" && rmm == 0 && tfb > 23 && tfb < 40 ) continue;
       if( det == "P0D" && rmm == 2 && tfb > 7 && tfb < 24 ) continue;

       for(int tript = 0; tript<=3; tript++){
   for(int chid = 0; chid<16; chid++){
     //	   std::cout << "ReferenceChannel: rmm = " << rmm << " tfb = " << tfb << "  tript = " << tript << "  channelID = " << chid << "  chan = " << chan << "  loped = " << ref_loped[rmm][tfb][tript][chid] <<std::endl;
     //	   std::cout << "ActualChannel: rmm = " << rmm << " tfb = " << tfb << "  tript = " << tript << "  channelID = " << chid << "  chan = " << chan << "  loped = " << loped[rmm][tfb][tript][chid] <<std::endl;

     chan = chid + tript*32 + tfb*4096 + rmm*524288 + firstCh;
     //	   std::cout << "Converted channel = " << chan << std::endl;

     //*****Skip channel if it is disconnected/unused
     bool unused_channel = false;
     for(unsigned int ch = 0; ch < NunusedCh; ch++){
       if( chan == unusedCh[ch] ){
         unused_channel = true;
         break;
       }
     }
     if( unused_channel ) continue;
     
     float lopedDrift = 0;
     float hipedDrift = 0;
     lopedDrift = loped[rmm][tfb][tript][chid] - ref_loped[rmm][tfb][tript][chid];
     hipedDrift = hiped[rmm][tfb][tript][chid] - ref_hiped[rmm][tfb][tript][chid];
     

     if (fabs(lopedDrift)>MaxDrift) {
       std::cout << "Timestamp: " << times[pedcounter] << " -- DRIFTING CHANNEL! rmm = " << rmm << "  tfb = " << tfb << "  tript = " << tript << "  channel = " << chid << "  Low Gain DRIFT " << lopedDrift << std::endl;
     }
     if (fabs(hipedDrift)>MaxDrift) {
       std::cout << "Timestamp: " << times[pedcounter] << " -- DRIFTING CHANNEL! rmm = " << rmm << "  tfb = " << tfb << "  tript = " << tript << "  channel = " << chid << "  High Gain DRIFT " << hipedDrift << std::endl;
     }

     //if chnannel drift is outside plot bounday
     if( hipedDrift > Yhigh ) hioverflow++;
     if( hipedDrift < Ylow ) hiunderflow++;
     if( lopedDrift > Yhigh ) looverflow++;
     if( lopedDrift < Ylow ) lounderflow++;
     
     
     if(loped[rmm][tfb][tript][chid]>MaxPed) {
       std::cout << "Timestamp: " << times[pedcounter] << " -- Low Gain channel with high pedestal value! rmm = " << rmm << "  tfb = " << tfb << "  tript = " << tript << "  channel = " << chid << "  Low Gain Ped =  " << loped[rmm][tfb][tript][chid] << std::endl;
     }
     if(loped[rmm][tfb][tript][chid]<MinPed) {
       std::cout << "Timestamp: " << times[pedcounter] << " -- Low Gain channel with low pedestal value! rmm = " << rmm << "  tfb = " << tfb << "  tript = " << tript << "  channel = " << chid << "  Low Gain Ped =  " << loped[rmm][tfb][tript][chid] << std::endl;
      }
     if(hiped[rmm][tfb][tript][chid]>MaxPed) {
       std::cout << "Timestamp: " << times[pedcounter] << " -- High Gain channel with high pedestal value! rmm = " << rmm << "  tfb = " << tfb << "  tript = " << tript << "  channel = " << chid << "  High Gain Ped =  " << hiped[rmm][tfb][tript][chid] << std::endl;
     }
     if(hiped[rmm][tfb][tript][chid]<MinPed) {
       std::cout << "Timestamp: " << times[pedcounter] << " --High Gain channel with low pedestal value! rmm = " << rmm << "  tfb = " << tfb << "  tript = " << tript << "  channel = " << chid << "  High Gain Ped =  " << hiped[rmm][tfb][tript][chid] << std::endl;
     }
     
     dptLoPedDriftPlot->Fill(times[pedcounter],((float)lopedDrift));
     dptHiPedDriftPlot->Fill(times[pedcounter],((float)hipedDrift));
     LoPedAll->Fill(times[pedcounter],((float)loped[rmm][tfb][tript][chid]));
     HiPedAll->Fill(times[pedcounter],((float)hiped[rmm][tfb][tript][chid]));

     dptLoPedDriftPlotRMM[rmm]->Fill(times[pedcounter],((float)lopedDrift));
     dptHiPedDriftPlotRMM[rmm]->Fill(times[pedcounter],((float)hipedDrift));
     dptLoPedPlotRMM[rmm]->Fill(times[pedcounter],((float)loped[rmm][tfb][tript][chid]));
     dptHiPedPlotRMM[rmm]->Fill(times[pedcounter],((float)hiped[rmm][tfb][tript][chid]));

          
   }// end loop over channel ID
       }// end loop over tript
     }// end loop over tfb
   } //end loop over rmm

   subtimes[pedcounter] = (int)(times[pedcounter]);
   hioverchsubrun[pedcounter] = hioverflow;
   hiunderchsubrun[pedcounter] = hiunderflow;
   looverchsubrun[pedcounter] = looverflow;
   lounderchsubrun[pedcounter] = lounderflow;

   pedcounter++;
       
 } //end of while(true)


 TGraph* HiOverChGraph = new TGraph(pedcounter,subtimes,hioverchsubrun);
 HiOverChGraph->SetName("HiOverChGraph");
 HiOverChGraph->SetTitle("High Overflow Channels");
 HiOverChGraph->SetLineColor(2);
 HiOverChGraph->Write();
 TGraph* HiUnderChGraph = new TGraph(pedcounter,subtimes,hiunderchsubrun);
 HiUnderChGraph->SetName("HiUnderChGraph");
 HiUnderChGraph->SetTitle("High Underflow Channels");
 HiUnderChGraph->SetLineColor(4);
 HiUnderChGraph->Write();
 TGraph* LoOverChGraph = new TGraph(pedcounter,subtimes,looverchsubrun);
 LoOverChGraph->SetName("LoOverChGraph");
 LoOverChGraph->SetTitle("Low Overflow Channels");
 LoOverChGraph->SetLineColor(2);
 LoOverChGraph->Write();
 TGraph* LoUnderChGraph = new TGraph(pedcounter,subtimes,lounderchsubrun);
 LoUnderChGraph->SetName("LoUnderChGraph");
 LoUnderChGraph->SetTitle("Low Underflow Channels");
 LoUnderChGraph->SetLineColor(4);
 LoUnderChGraph->Write();

 fOut.Write();

 //Plot and save histograms in ps file
 std::string psfile = Form("DQPedPlots%s.ps",det.c_str());
 std::string psfilefirst = psfile + "(";
 std::string psfilelast = psfile + ")";

 gStyle->SetOptStat(0);
 gStyle->SetPalette(1,0);

 // Pedestal for all RMMs combined
 TCanvas *c0 = new TCanvas("c0","",50,50,700,1000);
 c0->Divide(1,2);
 c0->cd(1);
 gPad->SetLogz();
 HiPedAll->SetMinimum(1);
 HiPedAll->Draw("colz");
 c0->cd(2);
 gPad->SetLogz();
 LoPedAll->SetMinimum(1);
 LoPedAll->Draw("colz");

 c0->Print(psfilefirst.c_str(),"Portrait");

 // Drift for all RMMs combined
 TCanvas *c0a = new TCanvas("c0","",50,50,700,1000);
 c0a->Divide(1,2);
 c0a->cd(1);
 gPad->SetLogz();
 dptHiPedDriftPlot->SetMinimum(1);
 dptHiPedDriftPlot->Draw("colz");
 c0a->cd(2);
 gPad->SetLogz();
 dptLoPedDriftPlot->SetMinimum(1);
 dptLoPedDriftPlot->Draw("colz");

 c0a->Print(psfile.c_str());

 // Over/underflow channel counts
 TCanvas *c0b = new TCanvas("c0b","",150,50,700,1000);
 c0b->Divide(1,2);
 c0b->cd(1);
 TMultiGraph *mg0 = new TMultiGraph();
 mg0->SetTitle("High gain: Num of over (red) and underflow (blue) chs");
 mg0->Add(HiOverChGraph,"lp");
 mg0->Add(HiUnderChGraph,"lp");
 mg0->Draw("a");
 mg0->GetXaxis()->SetTimeDisplay(1);
 mg0->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
 mg0->GetXaxis()->SetTimeOffset(0,"gmt");
 c0b->cd(2);
 TMultiGraph *mg1 = new TMultiGraph();
 mg1->SetTitle("Low gain: Num of over (red) and underflow (blue) chs");
 mg1->Add(LoOverChGraph,"lp");
 mg1->Add(LoUnderChGraph,"lp");
 mg1->Draw("a");
 mg1->GetXaxis()->SetTimeDisplay(1);
 mg1->GetXaxis()->SetTimeFormat("#splitline{%d/%m}{%H:%M}");
 mg1->GetXaxis()->SetTimeOffset(0,"gmt");

 c0b->Print(psfile.c_str());


 // Drift for individual RMMs: hi and lo per canvas
 TCanvas *crmm[nRMM];
 for( unsigned short int r = 0; r < nRMM; r++){
   std::string cname = Form("crmm%i",r);
   crmm[r] = new TCanvas(cname.c_str(),"",150,50,700,1000);
   crmm[r]->Divide(1,2);
   crmm[r]->cd(1);
   gPad->SetLogz();
   dptHiPedDriftPlotRMM[r]->SetMinimum(1);
   dptHiPedDriftPlotRMM[r]->Draw("colz");
   crmm[r]->cd(2);
   gPad->SetLogz();
   dptLoPedDriftPlotRMM[r]->SetMinimum(1);
   dptLoPedDriftPlotRMM[r]->Draw("colz");
   if( r < nRMM-1 ){
     crmm[r]->Print(psfile.c_str());
   } else {// last canvas
     crmm[r]->Print(psfilelast.c_str());
   }
 }

 return(0); 
}