source: trunk/FACT++/src/makeschedule.cc@ 18428

#include "externals/Prediction.h"

//#include <unordered_map>
#include <boost/algorithm/string/join.hpp>

#include "Database.h"

#include "tools.h"
#include "Time.h"
#include "Configuration.h"

/*
#include <TROOT.h>
#include <TH1.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TLegend.h>
*/

using namespace std;
using namespace Nova;

// -----------------------------------------------------------------------

/*
void CheckForGap(TCanvas &c, TGraph &g, double axis)
{
    if (g.GetN()==0 || axis-g.GetX()[g.GetN()-1]<450)
        return;

    c.cd();
    ((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
    while (g.GetN())
        g.RemovePoint(0);
}

void DrawClone(TCanvas &c, TGraph &g)
{
    if (g.GetN()==0)
        return;

    c.cd();
    ((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
}
*/

// ========================================================================
// ========================================================================
// ========================================================================

void SetupConfiguration(Configuration &conf)
{
    po::options_description control("Makeschedule");
    control.add_options()
        ("date", var<string>(), "SQL time (UTC), e.g. '2016-12-24'")
        ("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
        ("max-current", var<double>(90), "Global maximum current limit in uA")
        ("max-zd", var<double>(75), "Global zenith distance limit in degree")
        ("source", vars<string>(), "List of all TeV sources to be included, names according to the database")
        ("setup.*", var<double>(), "Setup for the sources to be observed")
        ("preobs.*", vars<string>(), "Prescheduled observations")
        ("startup.offset", var<double>(15), "Determines how many minutes the startup is scheduled before data-taking.start [0;120]")
        ("data-taking.start", var<double>(-12), "Begin of data-taking in degree of sun below horizon")
        ("data-taking.end", var<double>(-13.75), "End of data-taking in degree of sun below horizon")
        ;

    po::positional_options_description p;
    p.add("date", 1); // The first positional options

    conf.AddOptions(control);
    conf.SetArgumentPositions(p);
}

void PrintUsage()
{
    cout <<
        "makeschedule - Creates an automatic schedule\n"
        "\n"
        "Usage: makeschedule [yyyy-mm-dd]\n";
    cout << endl;
}

void PrintHelp()
{
#ifdef HAVE_ROOT
    cout <<
        "\n"
        "Examples:\n"
        "\n"
        "  makeschedule 2016-12-24\n"
        "\n"
        "Calculate the Christmas schedule for 2016 using all TeV sources from the\n"
        "database. If the date is omitted the current date is used.\n"
        "\n"
        "  makeschedule --source='Mrk 421' --source='Mrk 501' --source='Crab'\n"
        "\n"
        "Use only the mentioned sources to calculate the schedule.\n"
        "\n"
        "  makeschedule --source=Crab --setup.Crab.max-zd=30\n"
        "\n"
        "Limit the zenith distance of Crab into the range [0;30]deg.\n"
        "\n"
        "  makeschedule --source=Crab --setup.Crab.max-current=50\n"
        "\n"
        "Limit the maximum estimated current of Crab at 50uA.\n"
        "\n"
        "  makeschedule --source='IC 310' '--setup.IC 310.penalty=1.2'\n"
        "\n"
        "Multiply IC310's estimated relative threshold by a factor 1.2\n"
        "\n";
    cout << endl;
#endif
}


struct MyDouble
{
    double val;
    bool valid;
    MyDouble(Configuration &conf, const string &str) : val(0)
    {
        valid = conf.Has(str);
        if (valid)
            val = conf.Get<double>(str);
    }
    MyDouble() : val(0), valid(false) {}
};

/*
struct MinMax
{
    MyDouble min;
    MyDouble max;
    MinMax(Configuration &conf, const string &str)
    {
        min = MyDouble(conf, str+".min");
        max = MyDouble(conf, str+".max");
    }
    MinMax() {}
};
*/

struct Source
{
    // Global limits
    static double max_current;
    static double max_zd;

    // Source descrition
    string name;
    EquPosn equ;

    // Source specific limits
    MyDouble maxzd;
    MyDouble maxcurrent;
    double penalty;

    // Possible observation time
    double begin;
    double end;

    // Threshold (sorting reference)
    double threshold;

    double duration() const { return end-begin; };

    // Pre-observations (e.g. ratescan)
    vector<string> preobs;

    Source(const string &n="") : name(n), begin(0), threshold(std::numeric_limits<double>::max()) { }

    //bool IsSpecial() const { return threshold==std::numeric_limits<double>::max(); }

    double zd(const double &jd) const
    {
        return 90-GetHrzFromEqu(equ, jd).alt;
    }

    bool valid(const SolarObjects &so) const
    {
        const HrzPosn hrz     = GetHrzFromEqu(equ, so.fJD);
        const double  current = FACT::PredictI(so, equ);

        if (current>max_current)
            return false;

        if (hrz.alt<=0 || 90-hrz.alt>max_zd)
            return false;

        if (maxzd.valid && 90-hrz.alt>maxzd.val)
            return false;

        if (maxcurrent.valid && current>maxcurrent.val)
            return false;

        return true;
    }

    bool IsRangeValid(const double &jd_begin, const double &jd_end) const
    {
        const double step = 1./24/60;
        for (double jd=jd_begin; jd<jd_end+step/2; jd++)
            if (!valid(SolarObjects(jd)))
                return false;

        return true;
    }

    double getThreshold(const SolarObjects &so) const
    {
        const HrzPosn hrz = GetHrzFromEqu(equ, so.fJD);
        const double  current = FACT::PredictI(so, equ);

        const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);

        return penalty*ratio*pow(current/6.2, 0.394);
    }

    bool calcThreshold(const SolarObjects &so)
    {
        const HrzPosn hrz     = GetHrzFromEqu(equ, so.fJD);
        const double  current = FACT::PredictI(so, equ);

        if (current>max_current)
            return false;

        if (hrz.alt<=0 || 90-hrz.alt>max_zd)
            return false;

        if (maxzd.valid && 90-hrz.alt>maxzd.val)
            return false;

        if (maxcurrent.valid && current>maxcurrent.val)
            return false;

        const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
        threshold = penalty*ratio*pow(current/6.2, 0.394);

        return true;
    }
};

double Source::max_zd;
double Source::max_current;

bool SortByThreshold(const Source &i, const Source &j) { return i.threshold<j.threshold; }

bool RescheduleFirstSources(vector<Source> &obs)
{
    if (obs.size()<2 || obs[0].duration()>=40./24/60 || obs[0].name=="SLEEP" || obs[1].name=="SLEEP")
        return false;

    cout << "First source [" << obs[0].name << "] detected < 40min" << endl;

    const double obs1_duration = obs[1].end - obs[0].begin - 40./24/60;
    const double obs0_end      =              obs[0].begin + 40./24/60;

    // Check that:
    //  - the duration for the shrunken source obs[1] is still above 40min
    //  - obs[0] does not exceed 60deg at the end of its new window
    //  - obs[0] does not exceed any limit within the new window

    if (obs1_duration>=40./24/60 && obs[0].IsRangeValid(obs[0].end, obs0_end))
    {
        obs[0].end   = obs0_end;
        obs[1].begin = obs0_end;

        cout << "First source [" << obs[0].name << "] extended to 40min" << endl;

        return false;
    }

    // Try to remove the first source, check if second source fullfills all limits
    if (obs[1].IsRangeValid(obs[0].begin, obs[0].end))
    {
        cout << "First source [" << obs[0].name << "] removed" << endl;

        obs[1].begin = obs[0].begin;
        obs.erase(obs.begin());

        return true;
    }

    // Try to remove the second source, check if the first source fullfills all limits
    if (obs[0].IsRangeValid(obs[1].begin, obs[1].end))
    {
        cout << "Second source [" << obs[1].name << "] removed" << endl;

        obs[0].end = obs[1].end;
        obs.erase(obs.begin()+1);

        if (obs.size()==0 || obs[0].name!=obs[1].name)
            return true;

        obs[0].end = obs[1].end;
        obs.erase(obs.begin()+1);

        cout << "Combined first two indentical sources [" << obs[0].name << "] into one observation" << endl;

        return true;
    }

    cout << "No reschedule possible within limit." << endl;

    return false;
}

bool RescheduleLastSources(vector<Source> &obs)
{
    // If observation time is smaller than 40min for the first source
    // extend it to 40min if zenith angle will not go above 60deg.
    const int last = obs.size()-1;
    if (obs.size()<2 || obs[last].duration()>=40./24/60 || obs[last].name=="SLEEP" || obs[last-1].name=="SLEEP")
        return false;

    cout << "Last source [" << obs[last].name << "] detected < 40min" << endl;

    const double obs1_duration = obs[last].end - 40./24/60 - obs[last-1].begin;
    const double obs0_begin    = obs[last].end - 40./24/60;

    // Check that:
    //  - the duration for the shrunken source obs[1] is still above 40min
    //  - obs[0] does not exceed 60deg at the end of its new window
    //  - obs[0] does not exceed any limit within the new window

    if (obs1_duration>=40./24/60 && obs[last].IsRangeValid(obs0_begin, obs[last].begin))
    {
        obs[last].begin = obs0_begin;
        obs[last-1].end = obs0_begin;

        cout << "Last source [" << obs[last].name << "] extended to 40min" << endl;

        return false;
    }

    // Try to remove the last source, check if second source fullfills all limits
    if (obs[last-1].IsRangeValid(obs[last].begin, obs[last].end))
    {
        cout << "Last source [" << obs[last].name << "] removed" << endl;

        obs[last-1].end = obs[last].end;
        obs.pop_back();

        return true;
    }

    // Try to remove the second last source, check if the first source fullfills all limits
    if (obs[last].IsRangeValid(obs[last-1].begin, obs[last-1].end))
    {
        cout << "Second last source [" << obs[last-1].name << "] removed" << endl;

        obs[last].begin = obs[last-1].begin;
        obs.erase(obs.begin()+obs.size()-2);

        if (obs.size()==0 || obs[last-1].name!=obs[last-2].name)
            return true;

        obs[last-2].end = obs[last-1].end;
        obs.pop_back();

        cout << "Combined last two indentical sources [" << obs[last-1].name << "] into one observation" << endl;

        return true;
    }

    cout << "No reschedule possible within limit." << endl;

    return false;
}

bool RescheduleIntermediateSources(vector<Source> &obs)
{
    for (size_t i=1; i<obs.size()-1; i++)
    {
        if (obs[i].duration()>=40./24/60)
            continue;

        if (obs[i-1].name=="SLEEP" && obs[i+1].name=="SLEEP")
            continue;

        cout << "Intermediate source [" << obs[i].name << "] detected < 40min" << endl;

        double intersection = -1;

        if (obs[i-1].name=="SLEEP")
            intersection = obs[i].begin;

        if (obs[i+1].name=="SLEEP")
            intersection = obs[i].end;

        if (obs[i-1].name==obs[i+1].name)
            intersection = obs[i].begin;

        if (intersection<0)
        {
            const double step = 1./24/60;
            for (double jd=obs[i].begin; jd<obs[i].end-step/2; jd+=step)
            {
                const SolarObjects so(jd);
                if (obs[i-1].getThreshold(so)>=obs[i+1].getThreshold(so))
                {
                    intersection = jd;
                    break;
                }
            }
        }

        if ((obs[i-1].name!="SLEEP" && !obs[i-1].IsRangeValid(obs[i-1].end, intersection)) ||
            (obs[i+1].name!="SLEEP" && !obs[i+1].IsRangeValid(intersection, obs[i+1].begin)))
        {
            cout << "No reschedule possible within limits." << endl;
            continue;
        }

        cout << "Intermediate source [" << obs[i].name << "] removed" << endl;

        const bool underflow = obs[i-1].duration()*24*60<40 || obs[i+1].duration()*24*60<40;

        obs[i-1].end   = intersection;
        obs[i+1].begin = intersection;
        obs.erase(obs.begin()+i);

        i--;

        if (obs.size()>1 && obs[i].name==obs[i+1].name)
        {
            obs[i].end = obs[i+1].end;
            obs.erase(obs.begin()+i+1);

            cout << "Combined two surrounding indentical sources [" << obs[i].name << "] into one observation" << endl;

            i--;

            continue;
        }

        if (underflow)
            cout << "WARNING - Neighbor source < 40min as well." << endl;
    }
    return false;
}

void Print(const vector<Source> &obs, double startup_offset)
{
    cout << Time(obs[0].begin-startup_offset).GetAsStr() << "  STARTUP\n";
    for (const auto& src: obs)
    {
        string tm = Time(src.begin).GetAsStr();

        if (src.preobs.size()>0)
        {
            for (const auto& pre: src.preobs)
            {
                cout << tm << "  " << pre << "\n";
                tm = "                   ";
            }
        }

        cout << tm << "  " << src.name << " [";
        cout << src.duration()*24*60 << "'";
        if (src.name!="SLEEP")
            cout << Tools::Form("; %.1f/%.1f", src.zd(src.begin), src.zd(src.end));
        cout << "]";

        if (src.duration()*24*60<40)
            cout << " (!)";

        cout << "\n";
    }
    cout << Time(obs.back().end).GetAsStr() << "  SHUTDOWN" << endl;
}

int main(int argc, const char* argv[])
{
//    gROOT->SetBatch();

    Configuration conf(argv[0]);
    conf.SetPrintUsage(PrintUsage);
    SetupConfiguration(conf);

    if (!conf.DoParse(argc, argv, PrintHelp))
        return 127;

    // ------------------ Eval config ---------------------

    Time time;
    if (conf.Has("date"))
        time.SetFromStr(conf.Get<string>("date")+" 12:00:00");

    Source::max_current = conf.Get<double>("max-current");
    Source::max_zd      = conf.Get<double>("max-zd");

    const double startup_offset = conf.Get<double>("startup.offset")/60/24;

    const double angle_sun_set  = conf.Get<double>("data-taking.start");
    const double angle_sun_rise = conf.Get<double>("data-taking.end");

    if (startup_offset<0 || startup_offset>120)
        throw runtime_error("Only values [0;120] are allowed for startup.offset");

    if (angle_sun_set>-6)
        throw runtime_error("datataking.start not allowed before sun at -6deg");

    if (angle_sun_rise>-6)
        throw runtime_error("datataking.end not allowed after sun at -6deg");

    // -12: nautical
    // Sun set with the same date than th provided date
    // Sun rise on the following day
    const RstTime sun_set  = GetSolarRst(time.JD()-0.5, angle_sun_set);
    const RstTime sun_rise = GetSolarRst(time.JD()+0.5, angle_sun_rise);

    const double sunset  = ceil(sun_set.set*24*60)/24/60;
    const double sunrise = floor(sun_rise.rise*24*60)/24/60;

    cout << "\n";
    cout << "Date: " << Time(floor(sunset)) << "\n";
    cout << "Set:  " << Time(sun_set.set)   << "\n";
    cout << "Rise: " << Time(sun_rise.rise) << "\n";
    cout << endl;

    // ------------- Get Sources from databasse ---------------------

    const vector<string> ourcenames = conf.Vec<string>("source");
    const vector<string> sourcenames = conf.Vec<string>("source");
    cout << "Nsources = " << sourcenames.size() << "\n";

    string query = "SELECT fSourceName, fRightAscension, fDeclination FROM Source WHERE fSourceTypeKEY=1";
    if (sourcenames.size()>0)
        query += " AND fSourceName IN ('" + boost::algorithm::join(sourcenames, "', '")+"')";

    const string fDatabase = conf.Get<string>("source-database");
    const mysqlpp::StoreQueryResult res =
        Database(fDatabase).query(query).store();

    // ------------------ Eval config ---------------------

    vector<Source> sources;
    for (const auto &row: res)
    {
        const string name = string(row[0]);

        Source src(name);

        src.equ.ra  = double(row[1])*15;
        src.equ.dec = double(row[2]);

        src.maxzd = MyDouble(conf, "setup."+name+".max-zd");
        src.maxcurrent = MyDouble(conf, "setup."+name+".max-current");
        src.penalty = conf.Has("setup."+name+".penalty") ?
            conf.Get<double>("setup."+name+".penalty") : 1;

        src.preobs = conf.Vec<string>("preobs."+name);


        cout << "[" << name << "]";

        if (src.maxzd.valid)
            cout << " Zd<" << src.maxzd.val;
        if (src.penalty!=1)
            cout << " Penalty=" << src.penalty;

        cout << " " << boost::algorithm::join(src.preobs, "+") << endl;

        /*
         RstTime t1 = GetObjectRst(floor(sunset)-1, src.equ);
         RstTime t2 = GetObjectRst(floor(sunset),   src.equ);

         src.rst.transit = t1.transit<floor(sunset) ? t2.transit : t1.transit;
         src.rst.rise = t1.rise>src.rst.transit ? t2.rise : t1.rise;
         src.rst.set  = t1.set <src.rst.transit ? t2.set  : t1.set;
         */

        sources.emplace_back(src);
    }
    cout << endl;

    // -------------------------------------------------------------------------

    vector<Source> obs;

    const double step = 1./24/60;
    for (double jd=sunset; jd<sunrise-step/2; jd+=step)
    {
        const SolarObjects so(jd);

        vector<Source> vis;
        for (auto& src: sources)
        {
            if (src.calcThreshold(so))
                vis.emplace_back(src);
        }

        // In case no source was found, add a sleep source
        Source src("SLEEP");
        vis.emplace_back(src);

        // Source has higher priority if minimum observation time not yet fullfilled
        sort(vis.begin(), vis.end(), SortByThreshold);

        if (obs.size()>0 && obs.back().name==vis[0].name)
            continue;

        vis[0].begin = jd;
        obs.emplace_back(vis[0]);
    }

    if (obs.size()==0)
    {
        cout << "No source found." << endl;
        return 0;
    }

    // -------------------------------------------------------------------------

    for (auto it=obs.begin(); it<obs.end()-1; it++)
        it[0].end = it[1].begin;
    obs.back().end = sunrise;

    // -------------------------------------------------------------------------

    Print(obs, startup_offset);
    cout << endl;

    // -------------------------------------------------------------------------

    while (RescheduleFirstSources(obs));
    while (RescheduleLastSources(obs));
    while (RescheduleIntermediateSources(obs));

    // ---------------------------------------------------------------------

    cout << endl;
    Print(obs, startup_offset);
    cout << endl;

    // ---------------------------------------------------------------------

    return 1;

    // ------------- Create canvases and frames ---------------------
/*
    // It is important to use an offset which is larger than
    // 1970-01-01 00:00:00. This one will not work if your
    // local time zone is positive!
    TH1S hframe("", "", 1, Time(sunset).Mjd()*24*3600, Time(sunrise).Mjd()*24*3600);
    hframe.SetStats(kFALSE);
    hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
    hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+"  -  "+Time(jd+1).GetAsStr("%d/%m/%Y")+"  [UTC]").c_str());
    hframe.GetXaxis()->CenterTitle();
    hframe.GetYaxis()->CenterTitle();
    hframe.GetXaxis()->SetTimeDisplay(true);
    hframe.GetYaxis()->SetTitleSize(0.040);
    hframe.GetXaxis()->SetTitleSize(0.040);
    hframe.GetXaxis()->SetTitleOffset(1.1);
    hframe.GetYaxis()->SetLabelSize(0.040);
    hframe.GetXaxis()->SetLabelSize(0.040);

    TCanvas c1;
    c1.SetFillColor(kWhite);
    c1.SetBorderMode(0);
    c1.SetFrameBorderMode(0);
    c1.SetLeftMargin(0.085);
    c1.SetRightMargin(0.01);
    c1.SetTopMargin(0.03);
    c1.SetGrid();
    hframe.GetYaxis()->SetTitle("Altitude [deg]");
    hframe.SetMinimum(15);
    hframe.SetMaximum(90);
    hframe.DrawCopy();

    TCanvas c2;
    c2.SetFillColor(kWhite);
    c2.SetBorderMode(0);
    c2.SetFrameBorderMode(0);
    c2.SetLeftMargin(0.085);
    c2.SetRightMargin(0.01);
    c2.SetTopMargin(0.03);
    c2.SetGrid();
    hframe.GetYaxis()->SetTitle("Predicted Current [#muA]");
    hframe.SetMinimum(0);
    hframe.SetMaximum(100);
    hframe.DrawCopy();

    TCanvas c3;
    c3.SetFillColor(kWhite);
    c3.SetBorderMode(0);
    c3.SetFrameBorderMode(0);
    c3.SetLeftMargin(0.085);
    c3.SetRightMargin(0.01);
    c3.SetTopMargin(0.03);
    c3.SetGrid();
    c3.SetLogy();
    hframe.GetYaxis()->SetTitle("Estimated relative threshold");
    hframe.GetYaxis()->SetMoreLogLabels();
    hframe.SetMinimum(0.9);
    hframe.SetMaximum(11);
    hframe.DrawCopy();

    TCanvas c4;
    c4.SetFillColor(kWhite);
    c4.SetBorderMode(0);
    c4.SetFrameBorderMode(0);
    c4.SetLeftMargin(0.085);
    c4.SetRightMargin(0.01);
    c4.SetTopMargin(0.03);
    c4.SetGrid();
    hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
    hframe.SetMinimum(0);
    hframe.SetMaximum(180);
    hframe.DrawCopy();
    Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
    Int_t style[] = { kSolid, kDashed, kDotted };

    //    TLegend leg(0, 0, 1, 1);

    // ------------- Loop over sources ---------------------

    for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
    {
        // Eval row
        const string name = (*v)[0].c_str();

        EquPosn pos;
        pos.ra  = double((*v)[1])*15;
        pos.dec = double((*v)[2]);

        // Loop over 24 hours
        for (double h=0; h<1; h+=1./(24*12))
        {
            const SolarObjects so(jd+h);

            // get local position of source
            const HrzPosn hrz = GetHrzFromEqu(pos, so.fJD);

            if (v==res.begin())
                cout << Time(so.fJD) <<" " << 90-so.fMoonHrz.alt <<  endl;

            const double cur = FACT::PredictI(so, pos);

            // Relative  energy threshold prediction
            const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);

            // Add points to curve
            const double axis = Time(so.fJD).Mjd()*24*3600;

            // If there is a gap of more than one bin, start a new curve

            // Add data
            if (no_limits || cur<max_current)
                g1.SetPoint(g1.GetN(), axis, hrz.alt);

            if (no_limits || 90-hrz.alt<max_zd)
                g2.SetPoint(g2.GetN(), axis, cur);

            if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
                g3.SetPoint(g3.GetN(), axis, ratio*pow(cur/6.2, 0.394));

            if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
            {
                const double angle = GetAngularSeparation(so.fMoonEqu, pos);
                g4.SetPoint(g4.GetN(), axis, angle);
            }

            if (cnt==0)
                gm.SetPoint(gm.GetN(), axis, so.fMoonHrz.alt);
        }
    }

    // Save three plots
    TCanvas c5;
    c5.SetFillColor(kWhite);
    c5.SetBorderMode(0);
    c5.SetFrameBorderMode(0);
    leg.Draw();

    const string t = Time(jd).GetAsStr("%Y%m%d");

    c1.SaveAs((t+"-ZenithDistance.eps").c_str());
    c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
    c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
    c4.SaveAs((t+"-MoonDist.eps").c_str());
    c5.SaveAs((t+"-Legend.eps").c_str());

    c1.SaveAs((t+"-ZenithDistance.root").c_str());
    c2.SaveAs((t+"-PredictedCurrent.root").c_str());
    c3.SaveAs((t+"-RelativeThreshold.root").c_str());
    c4.SaveAs((t+"-MoonDist.root").c_str());

    c1.Print((t+".pdf(").c_str(), "pdf");
    c2.Print((t+".pdf" ).c_str(), "pdf");
    c3.Print((t+".pdf" ).c_str(), "pdf");
    c4.Print((t+".pdf" ).c_str(), "pdf");
    c5.Print((t+".pdf)").c_str(), "pdf");
*/
}