#define Nanoreader_cxx
#include "Nanoreader.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>


//   In a ROOT session, you can do:
//      Root > .L Nanoreader.C
//      Root > Nanoreader t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch

void Nanoreader::skyLoop()
{
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntries();
   Long64_t nbytes = 0, nb = 0;

   const double d2r = TMath::Pi()/180.;
   const double r2d = TMath::Pi()/180.;

   Double_t XRec[100000];
   Double_t YRec[100000];
   Double_t XPol[100000];
   Double_t YPol[100000];
   Double_t XmPol[100000];
   Double_t YmPol[100000];
   int iread = 0;
   int ireco = 0;
   int ilightning = 0;
   cout << "starting event loop" << endl;
   for (Long64_t jentry=0; jentry<nentries&&jentry<100000;jentry++) {
     //      Long64_t ientry = LoadTree(jentry);
     //      if (ientry < 0) break;
     nb = fChain->GetEntry(jentry);   nbytes += nb ;
     // test reconstrution (skip nan's)
     cout << "evt " << iread << ": " << reca_toutm << " " << reca_poutm << endl;
     if (reca_toutm != reca_toutm || reca_poutm != reca_poutm || 
	 reca_toutm<0.||reca_toutm>90.||reca_poutm<0.||reca_poutm>360.){
       cout << "entry " << iread << " badly reconstructed: " << reca_toutm << " " << reca_poutm << endl;
       cout << "                       " << reca_toutf << " " << reca_poutf << endl;
     } else {
       if(reca_toutm==0.||reca_poutm==0.||reca_toutm==90.) cout << "strange values " << reca_toutm << " " << reca_poutm << endl; 
       XRec[ireco]=360. + 180. - reca_poutf;
       if(XRec[ireco] > 360.) XRec[ireco] = XRec[ireco] - 360.;
       YRec[ireco]=90.-reca_toutf;
       XPol[ireco]=reca_toutf*cos(d2r*reca_poutf);
       YPol[ireco]=reca_toutf*sin(d2r*reca_poutf);
       XmPol[ireco]=reca_toutm*cos(d2r*reca_poutm);
       YmPol[ireco]=reca_toutm*sin(d2r*reca_poutm);
       //cout << "ireco,XPol,YPol " << ireco << " " << XPol[ireco] << " " << YPol[ireco] << endl;
       ireco++;
     }
     iread++;
   }
   cout << "Read " << iread << " evts from which " << ilightning << " presumed lightning stroke and " << ireco << " where properly reconstructed" << endl;
   const int nev2read = iread;
   const int nev2reco = ireco;

  // si j'ai bien compris, Joel mets le sud à 0 et tourne dans le sens trigo pour Phi...
  // J'ai l'habitude de lire une carte avec le nord en haut, le sud en bas l'est à droite et l'ouest à gauche (si, si!)
  // donc je change:
  //  TGraph *gPol = new TGraph(nev2read,XPol,YPol);
  Double_t XPolgeo[nev2reco];
  Double_t YPolgeo[nev2reco];
  Double_t XmPolgeo[nev2reco];
  Double_t YmPolgeo[nev2reco];
  for (int i=0;i<nev2reco;i++) {
    XPolgeo[i] =  YPol[i];
    YPolgeo[i] = -XPol[i];
    XmPolgeo[i] =  YmPol[i];
    YmPolgeo[i] = -XmPol[i];
  }

  TCanvas*  cSkyR = new TCanvas("cSkyR","Landscape",800,700);
  gStyle->SetOptStat(0);
  if (gDirectory->Get("hdum1")) {
     hdum1 = (TH2F*)gDirectory->Get("hdum1");
  } else {
    hdum1 = new TH2F("hdum1","Landscape",1,0.,360.,1,-5.,90.);
  }
  hdum1->Draw();
  TGraph *gRec = new TGraph(nev2reco,XRec,YRec);
  gRec->SetMarkerColor(2);
  gRec->SetMarkerStyle(24);
  // plot horizon line
  plotHorizonR(); 

  gRec->Draw("Psame");
  cSkyR->Print("cSkyR.pdf","pdf");


  TCanvas*  cSky = new TCanvas("cSky","Sky map",800,700);
  //gStyle->SetOptStat(0);
  if (gDirectory->Get("hdum2")) {
     hdum2 = (TH2F*)gDirectory->Get("hdum2");
  } else {
    hdum2 = new TH2F("hdum2","Sky map",1,-121.,121.,1,-121.,121.);
  }
  hdum2->Draw();
  TGraph *gPol = new TGraph(nev2reco,XPolgeo,YPolgeo);
  gPol->SetMarkerColor(2);
  gPol->SetMarkerStyle(24);
  // plot horizon line
  plotHorizon(); 
  // plot grid
  plotGrid();

  gPol->Draw("Psame");

  cSky->Print("cSky.pdf","pdf");

 
  TCanvas*  cSkym = new TCanvas("cSkym","Sky map in NOY coordinates",800,700);
  //gStyle->SetOptStat(0);
  hdum2->Draw();
  TGraph *gPolm = new TGraph(nev2reco,XmPolgeo,YmPolgeo);
  gPolm->SetMarkerColor(2);
  gPolm->SetMarkerStyle(24);
  
  // plot grid
  plotGrid();

  gPolm->Draw("Psame");

  cSkym->Print("cSkym.pdf","pdf");

}


void Nanoreader::timeLoop()
{
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntries();
   Long64_t nbytes = 0, nb = 0;
   int run_last = 0;
   double jd_last = 0;
   double now = 0;
   int run = 0;
   double timerange_low=99999999.;
   double timerange_end=0.;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb ;
      double now = dtt_djd;
      if(now<timerange_low)timerange_low=now;
      if(now>timerange_end)timerange_end=now;
   }
   int ndays = (int)(timerange_end-timerange_low+1);
   int nbins = ndays*24*6;
   cout << "Time scan over " << ndays << " days starting at djd=" << timerange_low << " and ending at djd=" << timerange_end << endl;
   TH1D *h1 = new TH1D("h1","Time of evts (per 10')",nbins,timerange_low,timerange_low+((double)ndays));
   TH1D *h2 = new TH1D("h2","Poisson",60,0.,60.);
   TH2D *h2d = new TH2D("h2d","log10(Delta t) vs time",nbins,timerange_low,timerange_low+((double)ndays),90,0.,3.);
   TH2D *htint = new TH2D("htint","integrated active time",nbins,timerange_low,timerange_low+((double)ndays),400,0.,400.);
   TH1D *h1druns = new TH1D("h1druns","runs vs time",nbins,timerange_low,timerange_low+((double)ndays));
   double integrated_time_days = 0.;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb ;
      double now = dtt_djd;
      //double now = tt_day*24.*3600.+tt_hour*3600.+tt_min*60.+tt_sec;
      //      cout << ientry << " ev=" << tt_ev << " now=" << now << " " << tt_day << " " << tt_hour << " " << tt_min << " " << tt_sec << endl;
      // if (Cut(ientry) < 0) continue;
      h1->Fill(now,1.);
      int tbin = h1druns->GetBin(now);
      h1druns->SetBinContent(tbin,(double)tt_run);
      if(jd_last > 0 && run_last == tt_run ) {
	double deltat = now - jd_last;
	//	cout << ientry << " ev=" << tt_ev << " jd=" << now << " deltat=" << deltat << endl;
	h2->Fill(deltat*24.*60.,1.);
	if(deltat > 0) {
	  h2d->Fill(now,log10(deltat*24.*60.),1.);
	  if(deltat < 20./(24.*60)) {
	    integrated_time_days += deltat;
	    htint->Fill(now,integrated_time_days,1.);
	  }
	}
      }      
      jd_last = now;
      run_last = tt_run;
   }

   TCanvas* cTimes = new TCanvas("cTimes","Time",10,10,1000,700);
   //gStyle->SetOptStat(0);
   cTimes->Divide(1,2);
   cTimes->cd(1);
   h1->DrawCopy();
   cTimes->cd(2);
   //   h2->Fit("expo","C","",0.,10.);
   h2->Fit("expo","L","",0.,15.);
   h2->DrawCopy("E");
   //cTimes_2->SetLogy(1);
   TCanvas* cIntTimes = new TCanvas("cIntTimes","Integrated active time",10,10,1000,700);
   gStyle->SetOptStat(0);
   htint->DrawCopy("box");
   char IT[80];
   sprintf(IT,"Total active time (in days) %7.3f",integrated_time_days);
   TLegend *legIT = new TLegend(0.1,0.85,0.9,0.9);
   legIT->SetHeader(IT);
   legIT->Draw();
   TCanvas* cTimes2d = new TCanvas("cTimes2d","Time 2D",100,100,1100,1100);
   cTimes2d->Divide(1,2);
   cTimes2d->cd(1);
   //gStyle->SetOptStat(0);
   cTimes2d->cd(1);
   h2d->DrawCopy("box");
   cTimes2d->cd(2);
   h1druns->Draw();   
}


void Nanoreader::tsLoop()
{
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();
   TH1D *hts = new TH1D("hts","ts adc bin",1024,0.,1024.);
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      int nb = fChain->GetEntry(jentry);
      double ts = d0_ts;
      if(d1_ts > ts) ts = d1_ts;
      if(d2_ts > ts) ts = d2_ts;
      if(d3_ts > ts) ts = d3_ts;
      if(d4_ts > ts) ts = d4_ts;
      ts=ts/4;
      cout << ientry << " " << ts << endl;
      hts->Fill(ts,1.);
   }

   TCanvas* cTS = new TCanvas("cTS","Trigger position",500,500);
   //gStyle->SetOptStat(0);
   //cTS->Divide(1,2);
   //cTS->cd(1);
   hts->Draw();
}



// to a set of RGB values.
void Nanoreader::Hue2RGB(float hue, float &r, float &g, float &b) {
  //assuming full saturation and mid intensity
  hue = hue + 360.*int(hue/360.);
  hue = hue - 360.*int(hue/360.);  
  //  float r,b,g;
  int n = int(hue/60.);
  const float m0 = 0.;
  const float m2 = 1.;
  float f = (hue -n*60)/60.;
  float mu = m0 +(m2-m0)*f;
  float md = m2 -(m2-m0)*f;
  //  cout << "n=" << n << endl;
  switch (n) {
  case 0 :	
    //Red-Yellow area
    r = m2;
    g = mu;
    b = m0;
    break;
  case 1 :	
    //Yellow-Green
    r = md;
    g = m2;
    b = m0;
    break;
  case 2 :  
    //Green-Cyan
    r = m0;
    g = m2;
    b = mu;
    break;
  case 3 :
    //Cyan-Blue
    r = m0;
    g = md;
    b = m2;
    break;
  case 4 :
    //Blue-Magenta
    r = mu;
    g = m0;
    b = m2;
    break;
  case 5 :
    //Magenta-Red
    r = m2;
    g = m0;
    b = md;
    break;
  default :
    r = 0;
    g = 0;
    b = 0;
  }
//   cout << "hue=" << hue;
//   cout << " r=" << r;
//   cout << " g=" << g;
//   cout << " b=" << b;
//   cout << endl;
}

void Nanoreader::plotQf(Long64_t entry, bool pdf)
{
  //0 1000 AHEAD30m50cm
  // NOY1-C 277633.00 5113127.00 511.9  1 1 1625.
  // NOY2   277645.27 5113099.08 512.0  2 2 1682.
  // NOY3   277660.99 5113139.12 512.0  3 3 1826.
  // NOY4   277621.29 5113155.16 512.0  4 4 1838.
  // NOY5   277604.67 5113115.71 512.0  5 5 1602.
  const int NSta = 5;
  const Double_t Easting[NSta] = {277633.00,277645.27,277660.99,277621.29,277604.67};
  const Double_t Northing[NSta] = {5113127.00,5113099.08,5113139.12,5113155.16,5113115.71};
  const Double_t Altitude[NSta] = {511.90,512.00,512.00,512.00,512.00};
  Double_t xet[NSta] = {0.,0.,0.,0.,0.};
  Double_t yet[NSta] = {0.,0.,0.,0.,0.};
  Double_t zet[NSta] = {0.,0.,0.,0.,0.};
  TCanvas* cHitMap;
  for(int is=0 ; is<NSta ; is++){
    xet[is]=(Easting[is]-Easting[0]) ;
    yet[is]=(Northing[is]-Northing[0]) ;
    zet[is]=Altitude[is]-Altitude[0] ;
  } // end is
  if (gDirectory->Get("hdummap")) {
    hdummap->Delete();
    //    cHitMap->Delete();
    gDetPos->Delete();
  }
  if (gDirectory->Get("cHitMap")) {
   //    cHitMap->Delete();
  } else {
    cHitMap = new TCanvas("cHitMap","Hit map",600,600);
  }
  gStyle->SetOptStat(0);
  Double_t w = gPad->GetWw()*gPad->GetAbsWNDC();
  Double_t h = gPad->GetWh()*gPad->GetAbsHNDC();
  Double_t xmin = 0;
  Double_t xmax = 2;
  Double_t ymin = 0;
  Double_t ymax = xmax*h/w; 
  cHitMap->SetFixedAspectRatio();
  cHitMap->Range(xmin,ymin,xmax,ymax); 
  gDetPos = new TGraph(NSta,xet,yet);
  gDetPos->SetMarkerColor(1);
  gDetPos->SetMarkerStyle(25);
  if (!fChain) return;
  fChain->GetEntry(entry);
  Double_t r0 = sqrt(d0_qf)*2.;
  Double_t r1 = sqrt(d1_qf)*2.;
  Double_t r2 = sqrt(d2_qf)*2.;
  Double_t r3 = sqrt(d3_qf)*2.;
  Double_t r4 = sqrt(d4_qf)*2.;
  Double_t rmax = 0.;
  if(r0>rmax) rmax=r0; if(r1>rmax) rmax=r1;  if(r2>rmax) rmax=r2;  if(r3>rmax) rmax=r3;  if(r4>rmax) rmax=r4;
  rmax=rmax+31;
  hdummap = new TH2F("hdummap","",1,-rmax,rmax,1,-rmax,rmax);
  hdummap->Draw();
  gDetPos->Draw("P");
  TEllipse *el0 = new TEllipse(xet[0],yet[0],r0,r0);
  TEllipse *el1 = new TEllipse(xet[1],yet[1],r1,r1);
  TEllipse *el2 = new TEllipse(xet[2],yet[2],r2,r2);
  TEllipse *el3 = new TEllipse(xet[3],yet[3],r3,r3);
  TEllipse *el4 = new TEllipse(xet[4],yet[4],r4,r4);
  Double_t trig = 999999.; Double_t tmax = 0.;
  Double_t ts0 = d0_ts; if (ts0<trig) trig=ts0; if (ts0>tmax) tmax=ts0;
  Double_t ts1 = d1_ts; if (ts1<trig) trig=ts1; if (ts1>tmax) tmax=ts1;
  Double_t ts2 = d2_ts; if (ts2<trig) trig=ts2; if (ts2>tmax) tmax=ts2;
  Double_t ts3 = d3_ts; if (ts3<trig) trig=ts3; if (ts3>tmax) tmax=ts3;
  Double_t ts4 = d4_ts; if (ts4<trig) trig=ts4; if (ts4>tmax) tmax=ts4;
  Double_t dtmax=tmax-trig;
  ts0=(ts0-trig)/dtmax;
  ts1=(ts1-trig)/dtmax;
  ts2=(ts2-trig)/dtmax;
  ts3=(ts3-trig)/dtmax;
  ts4=(ts4-trig)/dtmax;
  Float_t r_ = 0;
  Float_t g_ = 0;
  Float_t b_ = 0;
  // if (!gDirectory->Get("col0")) {
  //   col0->Delete();
  //   col1->Delete();
  //   col2->Delete();
  //   col3->Delete();
  //   col4->Delete();
  // }
  TColor *col0 = new TColor((Int_t)(20000+entry),0.,0.,0.);
  TColor *col1 = new TColor((Int_t)(21000+entry),0.,0.,0.);
  TColor *col2 = new TColor((Int_t)(22000+entry),0.,0.,0.);
  TColor *col3 = new TColor((Int_t)(23000+entry),0.,0.,0.);
  TColor *col4 = new TColor((Int_t)(24000+entry),0.,0.,0.);
  Float_t hue0 = ts0*300.;
  Float_t hue1 = ts1*300.;
  Float_t hue2 = ts2*300.;
  Float_t hue3 = ts3*300.;
  Float_t hue4 = ts4*300.;


  Hue2RGB(hue0,r_,g_,b_);cout << "d0 hue=" << hue0 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col0->SetRGB(r_,g_,b_);
  el0->SetFillColor(20000+entry);
  Hue2RGB(hue1,r_,g_,b_);cout << "d1 hue=" << hue1 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col1->SetRGB(r_,g_,b_);
  el1->SetFillColor(21000+entry);
  Hue2RGB(hue2,r_,g_,b_);cout << "d2 hue=" << hue2 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col2->SetRGB(r_,g_,b_);
  el2->SetFillColor(22000+entry);
  Hue2RGB(hue3,r_,g_,b_);cout << "d3 hue=" << hue3 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col3->SetRGB(r_,g_,b_);
  el3->SetFillColor(23000+entry);
  Hue2RGB(hue4,r_,g_,b_);cout << "d4 hue=" << hue4 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col4->SetRGB(r_,g_,b_);
  el4->SetFillColor(24000+entry);
  el0->Draw();
  el1->Draw();
  el2->Draw();
  el3->Draw();
  el4->Draw();
  //
  char text[80];
  TLatex ltext;
  ltext.SetTextAlign(12);
  ltext.SetTextSize(0.02);
  int n = 0;
  float xtxt = rmax/3;
  float ytxt = rmax/2;
  n = sprintf(text,"Entry = %lld",entry);
  ltext.DrawLatex(xtxt,1.7*ytxt,text);
  n = sprintf(text,"Run = %d Evt = %d",tt_run,tt_ev);
  ltext.DrawLatex(xtxt,1.6*ytxt,text);
  n = sprintf(text,"Date : %02d / %02d / %4d",tt_day,tt_month,tt_year);
  ltext.DrawLatex(xtxt,1.5*ytxt,text);
  n = sprintf(text,"Time : %02d:%02d:%02d",tt_hour,tt_min,tt_sec);
  ltext.DrawLatex(xtxt,1.4*ytxt,text);   
  if(reca_toutm == reca_toutm && reca_poutm == reca_poutm && 
     reca_toutm>0. && reca_toutm<90. && reca_poutm>0. && reca_poutm<360.)
    { 
      n = sprintf(text,"Reconstructed angles: ");
      ltext.DrawLatex(xtxt,1.3*ytxt,text);
      n = sprintf(text,"#theta_{m}=%5.2f  #phi_{m}=%5.2f",reca_toutm,reca_poutm);
      ltext.DrawLatex(xtxt,1.2*ytxt,text);
      n = sprintf(text,"#theta_{f}=%5.2f  #phi_{f}=%5.2f",reca_toutf,reca_poutf);
      ltext.DrawLatex(xtxt,1.1*ytxt,text);
    }
  n = sprintf(text,"d0.qf =%8.2e",d0_qf);
  ltext.DrawLatex(xet[0]+1,yet[0]+1,text);
  n = sprintf(text,"d1.qf =%8.2e",d1_qf);
  ltext.DrawLatex(xet[1]+1,yet[1]+1,text);
  n = sprintf(text,"d2.qf =%8.2e",d2_qf);
  ltext.DrawLatex(xet[2]+1,yet[2]+1,text);
  n = sprintf(text,"d3.qf =%8.2e",d3_qf);
  ltext.DrawLatex(xet[3]+1,yet[3]+1,text);
  n = sprintf(text,"d4.qf =%8.2e",d4_qf);
  ltext.DrawLatex(xet[4]+1,yet[4]+1,text);
  n = sprintf(text,"d0.ts =%8.2f ns",ts0*dtmax);
  ltext.DrawLatex(xet[0]+1,yet[0]-1,text);
  n = sprintf(text,"d1.ts =%8.2f ns",ts1*dtmax);
  ltext.DrawLatex(xet[1]+1,yet[1]-1,text);
  n = sprintf(text,"d2.ts =%8.2f ns",ts2*dtmax);
  ltext.DrawLatex(xet[2]+1,yet[2]-1,text);
  n = sprintf(text,"d3.ts =%8.2f ns",ts3*dtmax);
  ltext.DrawLatex(xet[3]+1,yet[3]-1,text);
  n = sprintf(text,"d4.ts =%8.2f ns",ts4*dtmax);
  ltext.DrawLatex(xet[4]+1,yet[4]-1,text);
  gDetPos->Draw("P");
  cHitMap->Update();
  if(pdf) {
    char title[80];
    int n = sprintf(title,"hitmap_%05lld.pdf",entry);
    cout << title << endl;
    cHitMap->Print(title,"pdf");
  }
  //
}

void Nanoreader::plotallQf()
{
  bool pdf = 0;
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   Long64_t nentries = fChain->GetEntries();
   for (Long64_t entry = 0; entry < nentries ; entry++) {
     plotQf(entry,pdf);
     // Wait for the return key to be press to continue
     pdf = 0;
     while (entry != nentries) {
       char c = getchar();
       if ( c == 'q'|| c == 'Q') entry = nentries;
       if ( c == 'p'|| c == 'P') {entry=entry-1; pdf=1; break;}
       if (c == '\n') break;

     }
   }
}

void Nanoreader::plotVEM(Long64_t entry, bool pdf)
{
  //0 1000 AHEAD30m50cm
  // NOY1-C 277633.00 5113127.00 511.9  1 1 1625.
  // NOY2   277645.27 5113099.08 512.0  2 2 1682.
  // NOY3   277660.99 5113139.12 512.0  3 3 1826.
  // NOY4   277621.29 5113155.16 512.0  4 4 1838.
  // NOY5   277604.67 5113115.71 512.0  5 5 1602.
  const int NSta = 5;
  const Double_t Easting[NSta] = {277633.00,277645.27,277660.99,277621.29,277604.67};
  const Double_t Northing[NSta] = {5113127.00,5113099.08,5113139.12,5113155.16,5113115.71};
  const Double_t Altitude[NSta] = {511.90,512.00,512.00,512.00,512.00};
  Double_t xet[NSta] = {0.,0.,0.,0.,0.};
  Double_t yet[NSta] = {0.,0.,0.,0.,0.};
  Double_t zet[NSta] = {0.,0.,0.,0.,0.};
  TCanvas* cHitMap;
  for(int is=0 ; is<NSta ; is++){
    xet[is]=(Easting[is]-Easting[0]) ;
    yet[is]=(Northing[is]-Northing[0]) ;
    zet[is]=Altitude[is]-Altitude[0] ;
  } // end is
  if (gDirectory->Get("hdummap")) {
    hdummap->Delete();
    //    cHitMap->Delete();
    gDetPos->Delete();
  }
  if (gDirectory->Get("cHitMap")) {
   //    cHitMap->Delete();
  } else {
    cHitMap = new TCanvas("cHitMap","Hit map",600,600);
  }
  gStyle->SetOptStat(0);
  Double_t w = gPad->GetWw()*gPad->GetAbsWNDC();
  Double_t h = gPad->GetWh()*gPad->GetAbsHNDC();
  Double_t xmin = 0;
  Double_t xmax = 2;
  Double_t ymin = 0;
  Double_t ymax = xmax*h/w; 
  cHitMap->SetFixedAspectRatio();
  cHitMap->Range(xmin,ymin,xmax,ymax); 
  gDetPos = new TGraph(NSta,xet,yet);
  gDetPos->SetMarkerColor(1);
  gDetPos->SetMarkerStyle(25);
  if (!fChain) return;
  fChain->GetEntry(entry);
  Double_t r0 = sqrt(d11_vem)*2.;
  Double_t r1 = sqrt(d21_vem)*2.;
  Double_t r2 = sqrt(d31_vem)*2.;
  Double_t r3 = sqrt(d41_vem)*2.;
  Double_t r4 = sqrt(d51_vem)*2.;
  Double_t rmax = 0.;
  if(r0>rmax) rmax=r0; if(r1>rmax) rmax=r1;  if(r2>rmax) rmax=r2;  if(r3>rmax) rmax=r3;  if(r4>rmax) rmax=r4;
  rmax=rmax+31;
  hdummap = new TH2F("hdummap","",1,-rmax,rmax,1,-rmax,rmax);
  hdummap->Draw();
  gDetPos->Draw("P");
  TEllipse *el0 = new TEllipse(xet[0],yet[0],r0,r0);
  TEllipse *el1 = new TEllipse(xet[1],yet[1],r1,r1);
  TEllipse *el2 = new TEllipse(xet[2],yet[2],r2,r2);
  TEllipse *el3 = new TEllipse(xet[3],yet[3],r3,r3);
  TEllipse *el4 = new TEllipse(xet[4],yet[4],r4,r4);
  Double_t trig = 999999.; Double_t tmax = 0.;
  Double_t ts0 = d0_ts; if (ts0<trig) trig=ts0; if (ts0>tmax) tmax=ts0;
  Double_t ts1 = d1_ts; if (ts1<trig) trig=ts1; if (ts1>tmax) tmax=ts1;
  Double_t ts2 = d2_ts; if (ts2<trig) trig=ts2; if (ts2>tmax) tmax=ts2;
  Double_t ts3 = d3_ts; if (ts3<trig) trig=ts3; if (ts3>tmax) tmax=ts3;
  Double_t ts4 = d4_ts; if (ts4<trig) trig=ts4; if (ts4>tmax) tmax=ts4;
  Double_t dtmax=tmax-trig;
  ts0=(ts0-trig)/dtmax;
  ts1=(ts1-trig)/dtmax;
  ts2=(ts2-trig)/dtmax;
  ts3=(ts3-trig)/dtmax;
  ts4=(ts4-trig)/dtmax;
  Float_t r_ = 0;
  Float_t g_ = 0;
  Float_t b_ = 0;
  // if (!gDirectory->Get("col0")) {
  //   col0->Delete();
  //   col1->Delete();
  //   col2->Delete();
  //   col3->Delete();
  //   col4->Delete();
  // }
  TColor *col0 = new TColor((Int_t)(20000+entry),0.,0.,0.);
  TColor *col1 = new TColor((Int_t)(21000+entry),0.,0.,0.);
  TColor *col2 = new TColor((Int_t)(22000+entry),0.,0.,0.);
  TColor *col3 = new TColor((Int_t)(23000+entry),0.,0.,0.);
  TColor *col4 = new TColor((Int_t)(24000+entry),0.,0.,0.);
  Float_t hue0 = ts0*300.;
  Float_t hue1 = ts1*300.;
  Float_t hue2 = ts2*300.;
  Float_t hue3 = ts3*300.;
  Float_t hue4 = ts4*300.;


  Hue2RGB(hue0,r_,g_,b_);cout << "d0 hue=" << hue0 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col0->SetRGB(r_,g_,b_);
  el0->SetFillColor(20000+entry);
  Hue2RGB(hue1,r_,g_,b_);cout << "d1 hue=" << hue1 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col1->SetRGB(r_,g_,b_);
  el1->SetFillColor(21000+entry);
  Hue2RGB(hue2,r_,g_,b_);cout << "d2 hue=" << hue2 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col2->SetRGB(r_,g_,b_);
  el2->SetFillColor(22000+entry);
  Hue2RGB(hue3,r_,g_,b_);cout << "d3 hue=" << hue3 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col3->SetRGB(r_,g_,b_);
  el3->SetFillColor(23000+entry);
  Hue2RGB(hue4,r_,g_,b_);cout << "d4 hue=" << hue4 << " rgb=" << r_ << "," << g_ << "," << b_ << endl;
  col4->SetRGB(r_,g_,b_);
  el4->SetFillColor(24000+entry);
  el0->Draw();
  el1->Draw();
  el2->Draw();
  el3->Draw();
  el4->Draw();
  //
  char text[80];
  TLatex ltext;
  ltext.SetTextAlign(12);
  ltext.SetTextSize(0.02);
  int n = 0;
  float xtxt = rmax/3;
  float ytxt = rmax/2;
  n = sprintf(text,"Entry = %lld",entry);
  ltext.DrawLatex(xtxt,1.7*ytxt,text);
  n = sprintf(text,"Run = %d Evt = %d",tt_run,tt_ev);
  ltext.DrawLatex(xtxt,1.6*ytxt,text);
  n = sprintf(text,"Date : %02d / %02d / %4d",tt_day,tt_month,tt_year);
  ltext.DrawLatex(xtxt,1.5*ytxt,text);
  n = sprintf(text,"Time : %02d:%02d:%02d",tt_hour,tt_min,tt_sec);
  ltext.DrawLatex(xtxt,1.4*ytxt,text);   
  if(reca_toutm == reca_toutm && reca_poutm == reca_poutm && 
     reca_toutm>0. && reca_toutm<90. && reca_poutm>0. && reca_poutm<360.)
    { 
      n = sprintf(text,"Reconstructed angles: ");
      ltext.DrawLatex(xtxt,1.3*ytxt,text);
      n = sprintf(text,"#theta_{m}=%5.2f  #phi_{m}=%5.2f",reca_toutm,reca_poutm);
      ltext.DrawLatex(xtxt,1.2*ytxt,text);
      n = sprintf(text,"#theta_{f}=%5.2f  #phi_{f}=%5.2f",reca_toutf,reca_poutf);
      ltext.DrawLatex(xtxt,1.1*ytxt,text);
    }
  n = sprintf(text,"d11.vem =%8.2e",d11_vem);
  ltext.DrawLatex(xet[0]+1,yet[0]+1,text);
  n = sprintf(text,"d21.vem =%8.2e",d21_vem);
  ltext.DrawLatex(xet[1]+1,yet[1]+1,text);
  n = sprintf(text,"d31.vem =%8.2e",d31_vem);
  ltext.DrawLatex(xet[2]+1,yet[2]+1,text);
  n = sprintf(text,"d41.vem =%8.2e",d41_vem);
  ltext.DrawLatex(xet[3]+1,yet[3]+1,text);
  n = sprintf(text,"d51.vem =%8.2e",d51_vem);
  ltext.DrawLatex(xet[4]+1,yet[4]+1,text);
  n = sprintf(text,"d0.ts =%8.2f ns",ts0*dtmax);
  ltext.DrawLatex(xet[0]+1,yet[0]-1,text);
  n = sprintf(text,"d1.ts =%8.2f ns",ts1*dtmax);
  ltext.DrawLatex(xet[1]+1,yet[1]-1,text);
  n = sprintf(text,"d2.ts =%8.2f ns",ts2*dtmax);
  ltext.DrawLatex(xet[2]+1,yet[2]-1,text);
  n = sprintf(text,"d3.ts =%8.2f ns",ts3*dtmax);
  ltext.DrawLatex(xet[3]+1,yet[3]-1,text);
  n = sprintf(text,"d4.ts =%8.2f ns",ts4*dtmax);
  ltext.DrawLatex(xet[4]+1,yet[4]-1,text);
  gDetPos->Draw("P");
  cHitMap->Update();
  if(pdf) {
    char title[80];
    int n = sprintf(title,"hitmap_%05lld.pdf",entry);
    cout << title << endl;
    cHitMap->Print(title,"pdf");
  }
  //
}

void Nanoreader::plotallVEM()
{
  bool pdf = 0;
// Print contents of entry.
// If entry is not specified, print current entry
   if (!fChain) return;
   Long64_t nentries = fChain->GetEntries();
   for (Long64_t entry = 0; entry < nentries ; entry++) {
     plotVEM(entry,pdf);
     // Wait for the return key to be press to continue
     pdf = 0;
     while (entry != nentries) {
       char c = getchar();
       if ( c == 'q'|| c == 'Q') entry = nentries;
       if ( c == 'p'|| c == 'P') {entry=entry-1; pdf=1; break;}
       if (c == '\n') break;

     }
   }
}

void Nanoreader::makeKml(Long64_t entry)
{
  const float deg2rad = 3.141592653589793/180.;
  const float lonAC = 6.121161;
  const float latAC = 46.135455;
  const float lonA1 = 6.121330;
  const float latA1 = 46.135189;
  const float lonA2 = 6.121530;
  const float latA2 = 46.135572;
  const float lonA3 = 6.120979;
  const float latA3 = 46.135711;
  const float lonA4 = 6.120779;
  const float latA4 = 46.135337;
  const float lon0 = lonAC;
  const float lat0 = latAC;

   if (!fChain) return;
   fChain->GetEntry(entry);
   bool reco = 1;
   bool lightning = 0;
   Double_t XRec[1];
   Double_t YRec[1];
   if(reca_toutm != reca_toutm || reca_poutm != reca_poutm || 
     reca_toutm<0.||reca_toutm>90.||reca_poutm<0.||reca_poutm>360.){
     cout << "evt " << entry << " badly reconstructed: " << reca_toutm << " " << reca_poutm << endl;
     reco = 0;
     lightning = 0;
   }
//

   if(reco == 1 && lightning == 0) {
     float theta = reca_toutf;
     float phi = 360. + 180. - reca_poutf;
     if(phi > 360.) phi = phi - 360.;
     cout << "phi (azimutal angle or heading) = " << phi << "  theta (zenith angle or tilt) = " << theta << endl;
     char kmlfile[80];
     int n = sprintf(kmlfile,"EvtDisp_%05lld.kml",entry);
     ofstream outfile;
     outfile.open(kmlfile);
     outfile << "<?xml version=\"1.0\" encoding=\"UTF-8\"?> <kml xmlns=\"http://www.opengis.net/kml/2.2\" xmlns:gx=\"http://www.google.com/kml/ext/2.2\" xmlns:kml=\"http://www.opengis.net/kml/2.2\" xmlns:atom=\"http://www.w3.org/2005/Atom\">";
     outfile << endl;
     outfile << "<Document>" << endl;
     outfile << "	<open>0</open>" << endl;
     outfile << "       <visibility>1</visibility>" << endl;
     outfile << "	<name>AHEAD Event display</name>" << endl;
     outfile << "	<description>Created by FM</description>" << endl;
     outfile << "	<LookAt id=\"camera_pos\">" << endl;
     outfile << "		<longitude>" << setprecision(10) << lon0 << "</longitude>" << endl;
     outfile << "		<latitude>"  << setprecision(10) << lat0 << "</latitude>" << endl;
     outfile << "		<altitude>0</altitude>" << endl;
     float camera_heading = phi+90.;
     if(camera_heading>360.) camera_heading = camera_heading - 360.;
     outfile << "		<heading>" << setprecision(10) << camera_heading << "</heading>" << endl;
     outfile << "		<tilt>70.</tilt>" << endl;
     outfile << "		<range>300.</range>" << endl;
     outfile << "	</LookAt>" << endl;
     outfile << "	<Style id=\"s_hit\">" << endl;
     outfile << "		<IconStyle id=\"is_hit\">" << endl;
     outfile << "			<scale>2</scale>" << endl;
     outfile << "			<color>ffffff55</color>" << endl;
     outfile << "			<Icon>" << endl;
     outfile << "				<href>hit.png</href>" << endl;
     outfile << "			</Icon>" << endl;
     outfile << "			<hotSpot x=\"0.5\" y=\"0.5\" xunits=\"fraction\" yunits=\"fraction\"/>" << endl;
     outfile << "		</IconStyle>" << endl;
     outfile << "		<ListStyle>" << endl;
     outfile << "		</ListStyle>" << endl;
     outfile << "	</Style>" << endl;
     const float dstep =  1; // 1m steps (3ns)
     const float tstart = 0;
     const float dstart = 50; // starting animation 100m away from AC
     const float t0offset = dstart*3; // time offset
     const float dstop = -30; // Ending underground 30m away from AC
     const float tstop = tstart+dstop*3+t0offset;
     const float tendanim = tstop+200; // Ending time for animation
     outfile << "	<Placemark>" << endl;
     outfile << "		<TimeSpan>" << endl;
     outfile << "			<begin>" << t0offset << "</begin>" << endl;
     outfile << "			<end>" << tendanim << "</end>" << endl;
     outfile << "		</TimeSpan>" << endl;
     outfile << "		<styleUrl>#s_hit</styleUrl>" << endl;
     outfile << "		<Point>" << endl;
     outfile << "			<coordinates>" << setprecision(10) << lonAC << "," << latAC << ",0</coordinates>" << endl;
     outfile << "		</Point>" << endl;
     outfile << "	</Placemark>" << endl;
     outfile << "	<Placemark>" << endl;
     outfile << "		<TimeSpan>" << endl;
     outfile << "			<begin>" << t0offset + d1_ts - d0_ts << "</begin>" << endl;
     outfile << "			<end>" << tendanim << "</end>" << endl;
     outfile << "		</TimeSpan>" << endl;
     outfile << "		<styleUrl>#s_hit</styleUrl>" << endl;
     outfile << "		<Point>" << endl;
     outfile << "			<coordinates>" << setprecision(10) << lonA1 << "," << latA1 << ",0</coordinates>" << endl;
     outfile << "		</Point>" << endl;
     outfile << "	</Placemark>" << endl;
     outfile << "	<Placemark>" << endl;
     outfile << "		<TimeSpan>" << endl;
     outfile << "			<begin>" << t0offset + d2_ts - d0_ts << "</begin>" << endl;
     outfile << "			<end>" << tendanim << "</end>" << endl;
     outfile << "		</TimeSpan>" << endl;
     outfile << "		<styleUrl>#s_hit</styleUrl>" << endl;
     outfile << "		<Point>" << endl;
     outfile << "			<coordinates>" << setprecision(10) << lonA2 << "," << latA2 << ",0</coordinates>" << endl;
     outfile << "		</Point>" << endl;
     outfile << "	</Placemark>" << endl;
     outfile << "	<Placemark>" << endl;
     outfile << "		<TimeSpan>" << endl;
     outfile << "			<begin>" << t0offset + d3_ts - d0_ts << "</begin>" << endl;
     outfile << "			<end>" << tendanim << "</end>" << endl;
     outfile << "		</TimeSpan>" << endl;
     outfile << "		<styleUrl>#s_hit</styleUrl>" << endl;
     outfile << "		<Point>" << endl;
     outfile << "			<coordinates>" << setprecision(10) << lonA3 << "," << latA3 << ",0</coordinates>" << endl;
     outfile << "		</Point>" << endl;
     outfile << "	</Placemark>" << endl;
     outfile << "	<Placemark>" << endl;
     outfile << "		<TimeSpan>" << endl;
     outfile << "			<begin>" << setprecision(4) << t0offset + d4_ts - d0_ts << "</begin>" << endl;
     outfile << "			<end>" << setprecision(4) << tendanim << "</end>" << endl;
     outfile << "		</TimeSpan>" << endl;
     outfile << "		<styleUrl>#s_hit</styleUrl>" << endl;
     outfile << "		<Point>" << endl;
     outfile << "			<coordinates>" << setprecision(10) << lonA4 << "," << latA4 << ",0</coordinates>" << endl;
     outfile << "		</Point>" << endl;
     outfile << "	</Placemark>" << endl;
     for (float dist = dstart; dist > dstop ; dist = dist - dstep) {
       outfile << "	<Placemark>" << endl;
       outfile << "		<name>sh_front</name>" << endl;
       outfile << "		<TimeSpan>" << endl;
       outfile << "			<begin>" << setprecision(4) << t0offset - dist*3 << "</begin>" << endl;
       outfile << "			<end>" <<  setprecision(4) << t0offset - dist*3 + 3.*dstep << "</end>" << endl;
       outfile << "		</TimeSpan>" << endl;
       outfile << "		<Model>" << endl;
       outfile << "			<altitudeMode>relativeToGround</altitudeMode>" << endl;
       const float alt = 0+dist*cos(theta*deg2rad);
       const float dNord = dist*sin(theta*deg2rad)*cos(phi*deg2rad);
       const float dEast = dist*sin(theta*deg2rad)*sin(phi*deg2rad);
       const float lat =  lat0 + dNord * (360./4.e7);
       const float lon =  lon0 + dEast * (360./4.e7)*cos(lat*deg2rad);
       outfile << "			<Location>" << endl;
       outfile << "		                <latitude>"  << setprecision(10) << lat << "</latitude>" << endl;
       outfile << "				<longitude>" << setprecision(10) << lon << "</longitude>" << endl;
       outfile << "				<altitude>"  << alt << "</altitude>" << endl;
       outfile << "			</Location>" << endl;
       outfile << "			<Orientation>" << endl;
       outfile << "				<heading>" << setprecision(8) << phi << "</heading>" << endl;
       outfile << "				<tilt>" << setprecision(8) << theta << "</tilt>" << endl;
       outfile << "       				<roll>0</roll>" << endl;
       outfile << "				</Orientation>" << endl;
       outfile << "			<Scale>" << endl;
       outfile << "				<x>1</x>" << endl;
       outfile << "				<y>1</y>" << endl;
       outfile << "				<z>1</z>" << endl;
       outfile << "			</Scale>" << endl;
       outfile << "			<Link>" << endl;
       outfile << "				<href>pancake.dae</href>" << endl;
       outfile << "			</Link>" << endl;
       outfile << "			<ResourceMap>" << endl;
       outfile << "			</ResourceMap>" << endl;
       outfile << "		</Model>" << endl;
       outfile << "	</Placemark>" << endl;
     }
     outfile << "</Document>" << endl;
     outfile << "</kml>" << endl;
   // close the opened file.
     outfile.close();
   }
}

void Nanoreader::makeJson(Long64_t entry)
{

   if (!fChain) return;
   fChain->GetEntry(entry);
   bool reco = 1;
   bool lightning = 0;
   if(reca_toutm != reca_toutm || reca_poutm != reca_poutm || 
     reca_toutm<0.||reca_toutm>90.||reca_poutm<0.||reca_poutm>360.){
     cout << "evt " << entry << " badly reconstructed: " << reca_toutm << " " << reca_poutm << endl;
     reco = 0;
   }
//

   if(reco == 1 && lightning == 0) {
     float theta = reca_toutf;
     float phi = 360. + 180. - reca_poutf;
     if(phi > 360.) phi = phi - 360.;
     cout << "phi (azimutal angle or heading) = " << phi << "  theta (zenith angle or tilt) = " << theta << endl;
     char Jsonfile[80];
     int n = sprintf(Jsonfile,"Json/AHEAD_Evt_%04d_%05lld.Json",tt_run,entry);
     ofstream outfile;
     outfile.open(Jsonfile);
     char date[80];
     n = sprintf(date,"%04d-%02d-%02dT%02d:%02d:%02d.%03dZ",tt_year,tt_month,tt_day,tt_hour,tt_min,tt_sec,0);
     cout << "{" << endl;
     cout << "   \"RunId\" : [" << tt_run << "]," << endl;
     cout << "   \"EvtId\" : [" << tt_ev  << "]," << endl;
     cout << "   \"EvtJD\" : \"" << date << "\"," << endl;
     cout << "   \"DetHitTimes\" : [" << d0_tf << "," << d1_tf << "," << d2_tf << "," << d3_tf << "," << d4_tf << "]," << endl;
     cout << "   \"xCore\" : [" << ldf_core[0] << "]," << endl;
     cout << "   \"yCore\" : [" << ldf_core[1] << "]," << endl;
     cout << "   \"Theta\" : [" << theta << "]," << endl;
     cout << "   \"Phi\" : [" << phi << "]" << endl;
     cout << "   \"ShSize\" : [" << ldf_Ne << "]" << endl;
     cout << "   \"Energy\" : [" << ldf_Ep << "]" << endl;
     cout << "}" << endl;
     //
     outfile << "{" << endl;
     outfile << "   \"RunId\" : [" << tt_run << "]," << endl;
     outfile << "   \"EvtId\" : [" << tt_ev  << "]," << endl;
     outfile << "   \"EvtJD\" : \"" << date << "\"," << endl;
     outfile << "   \"DetHitTimes\" : [" << d0_tf << "," << d1_tf << "," << d2_tf << "," << d3_tf << "," << d4_tf << "]," << endl;
     outfile << "   \"xCore\" : [" << ldf_core[0] << "]," << endl;
     outfile << "   \"yCore\" : [" << ldf_core[1] << "]," << endl;
     outfile << "   \"Theta\" : [" << theta << "]," << endl;
     outfile << "   \"Phi\" : [" << phi << "]," << endl;
     outfile << "   \"ShSize\" : [" << ldf_Ne << "]," << endl;
     outfile << "   \"Energy\" : [" << ldf_Ep << "]" << endl;
     outfile << "}" << endl;
   // close the opened file.
     outfile.close();
   }
}



void Nanoreader::plotCIC(double lambdaAtt){
  Long64_t nentries = fChain->GetEntriesFast();
  TH1D * h1 = new TH1D("h1","CIC 1.0<sec1.2",50,4.,9.);
  TH1D * h2 = new TH1D("h2","CIC 1.2<sec1.4",50,4.,9.);
  TH1D * h3 = new TH1D("h3","CIC 1.4<sec1.6",50,4.,9.);
  TH1D * h4 = new TH1D("h4","CIC 1.6<sec1.8",50,4.,9.);
  TH1D * h5 = new TH1D("h5","CIC 1.8<sec2.0",50,4.,9.);
  TH1D * hCIC1 = new TH1D("hCIC1","CIC 1.0<sec1.2",50,4.,9.);
  TH1D * hCIC2 = new TH1D("hCIC2","CIC 1.2<sec1.4",50,4.,9.);
  TH1D * hCIC3 = new TH1D("hCIC3","CIC 1.4<sec1.6",50,4.,9.);
  TH1D * hCIC4 = new TH1D("hCIC4","CIC 1.6<sec1.8",50,4.,9.);
  TH1D * hCIC5 = new TH1D("hCIC5","CIC 1.8<sec2.0",50,4.,9.);
  hCIC1->SetLineColor(1);
  hCIC2->SetLineColor(2);
  hCIC3->SetLineColor(3);
  hCIC4->SetLineColor(4);
  hCIC5->SetLineColor(6);
  
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    int nb = fChain->GetEntry(jentry);
    cout << "entry=" << jentry << endl;
    double sectheta = 1./cos(3.141592654*reca_toutf/180.);
    double log10Ne = log10(ldf_Ne);
    const double Xvert = 975.;  // atm depth for vertical at 500m
    //    double lambdaAtt 90. = .; // attenuation at the Nancay site
    // selection
    double rcore = sqrt(ldf_core[0]*ldf_core[0]+ldf_core[1]*ldf_core[1]);
    //    if(rcore < 35.) {
    if(reca_toutf<50.&&rcore < 35.&&ldf_Ecore[0]<10&&ldf_Ecore[1]<10) {
      double corlog10Ne = log10(ldf_Ne*exp(Xvert*(sectheta - 1.)/lambdaAtt));
      if(sectheta>=1.0&&sectheta<1.1) {
	h1->Fill(corlog10Ne);
      }else if(sectheta<1.2) {
	h2->Fill(corlog10Ne);
      }else if(sectheta<1.3) {
	h3->Fill(corlog10Ne);
      }else if(sectheta<1.4) {
	h4->Fill(corlog10Ne);
      }else if(sectheta<1.5) {
	h5->Fill(corlog10Ne);
      }
    }
  }
  // make integral spectra
  for(int i = 0; i<50; i++) {
    double integ = h1->GetBinContent(i);
    for (int j = i; j<50; j++) {
      integ += h1->GetBinContent(j);
    }
    hCIC1->SetBinContent(i,integ);
  }
  for(int i = 0; i<50; i++) {
    double integ = h2->GetBinContent(i);
    for (int j = i; j<50; j++) {
      integ += h2->GetBinContent(j);
    }
    hCIC2->SetBinContent(i,integ);
  }
  for(int i = 0; i<50; i++) {
    double integ = h3->GetBinContent(i);
    for (int j = i; j<50; j++) {
      integ += h3->GetBinContent(j);
    }
    hCIC3->SetBinContent(i,integ);
  }
  for(int i = 0; i<50; i++) {
    double integ = h4->GetBinContent(i);
    for (int j = i; j<50; j++) {
      integ += h4->GetBinContent(j);
    }
    hCIC4->SetBinContent(i,integ);
  }
  for(int i = 0; i<50; i++) {
    double integ = h5->GetBinContent(i);
    for (int j = i; j<50; j++) {
      integ += h5->GetBinContent(j);
    }
    hCIC5->SetBinContent(i,integ);
  }
   TCanvas* cCIC = new TCanvas("cCIC","CIC plots",500,500);
   gStyle->SetOptStat(0);
   //cCIC->Divide(1,2);
   //cCIC->cd(1);
   gPad->SetLogy(1);
   hCIC1->Draw();
   hCIC1->GetXaxis()->SetTitle("N*_{e}");
   hCIC1->GetYaxis()->SetTitle("Intensity(>N*_{e})");
   hCIC1->SetTitle("");
   hCIC2->Draw("same");
   hCIC3->Draw("same");
   hCIC4->Draw("same");
   hCIC5->Draw("same");
   TLegend *legCIC = new TLegend(0.6,0.5,0.99,0.99);
   legCIC->SetHeader("CIC method : intensity(shower_size > N*_{e}) vs N*_{e}");
   legCIC->AddEntry((TObject*)0,"N*_{e} = N_{e} #times exp(X_{vert}*(sec(#theta) - 1)/#lambda_{Att})","");
   char sl[80];
   sprintf(sl,"with #lambda_{Att} = %4.1f g/cm^2",lambdaAtt);
   legCIC->AddEntry((TObject*)0,sl,"");
   legCIC->AddEntry((TObject*)0,"and N_{e} from LDF at core distance = 0","");
   
   legCIC->AddEntry(hCIC1,"1.0 #leq sec(#theta) <1.1","l");
   legCIC->AddEntry(hCIC2,"1.1 #leq sec(#theta) <1.2","l");
   legCIC->AddEntry(hCIC3,"1.2 #leq sec(#theta) <1.3","l");
   legCIC->AddEntry(hCIC4,"1.3 #leq sec(#theta) <1.4","l");
   legCIC->AddEntry(hCIC5,"1.4 #leq sec(#theta) <1.5","l");
   legCIC->Draw();
}

void Nanoreader::plotSpectrum(double lambdaAtt){
  const double texp = 183.225*24.*3600.; //seconds
  const double area = 2.*pow(30.5,2.); //m²
  const double strad =  3.141592654; // down to 60°
  const double norma = 1./texp/area/strad/(pow(10.,0.1)-1.);
    
  Long64_t nentries = fChain->GetEntriesFast();
  TH1D * hspectrum = new TH1D("hspectrum","spectrum",33,15.,18.3);
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    int nb = fChain->GetEntry(jentry);
    cout << "entry=" << jentry << endl;
    double rcore = sqrt(ldf_core[0]*ldf_core[0]+ldf_core[1]*ldf_core[1]);
    if(reca_toutf<60.&&rcore < 30.5&&ldf_Ecore[0]<10&&ldf_Ecore[1]<10) {
      double sectheta = 1./cos(3.141592654*reca_toutf/180.);
      double log10Ne = log10(ldf_Ne);
      const double Xvert = 975.;  // atm depth for vertical at 500m
      //    double lambdaAtt = 90.; // attenuation at the Nancay site
      double corlog10Ne = log10(ldf_Ne*exp(Xvert*(sectheta - 1.)/lambdaAtt));
      double log10ECIC =  log10(2.183E10) + 0.9*corlog10Ne;
      double ECIC2 = pow(10.,2.*log10ECIC);
      double ECIC21 = pow(10.,2.1*log10ECIC);
      double ECIC = pow(10.,log10ECIC);
      double ECIC15 = pow(10.,1.5*log10ECIC);
      hspectrum->Fill(log10ECIC,norma*ECIC15);
    }
  }
  TCanvas* cECIC = new TCanvas("cECIC","ECIC spectrum",500,500);
  gStyle->SetOptStat(0);
  gPad->SetLogy(1);
  hspectrum->Draw();
  hspectrum->GetXaxis()->SetTitle("log10(ECIC)");
  hspectrum->GetYaxis()->SetTitle("E^{2.5}.J/(m^{-2}.s^{-1}.sr^{-1}.eV^{1.5})");
  hspectrum->SetMaximum(1.E17);
  hspectrum->SetMinimum(2.E14);
  hspectrum->SetLineColor(kMagenta);
  hspectrum->SetLineWidth(3);
}


void Nanoreader::plotLDF(Long64_t entry, bool pdf) {
   if (!fChain) return;
   fChain->GetEntry(entry);
   double Rmax = 0;
   float Er[5];
   float EQ[5];
   for (int is=0;  is<5; is++) {
     if(ldf_r[is] > Rmax) Rmax = ldf_r[is];
     Er[is] = sqrt(ldf_Ecore[0]*ldf_Ecore[0]+ldf_Ecore[1]*ldf_Ecore[1]);
     EQ[is] = sqrt(ldf_Q[is])*1.;
   }

   TF1 *fNKG = new TF1("fNKG","([0]*0.443/(6400.))*pow(x/80.,-0.8)*pow(1.+x/80.,-2.5)",0.,Rmax*1.1);

  fNKG->SetParameter(0,ldf_Ne);
  TCanvas* cLDF = new TCanvas("cLDF","LDF fit",500,500);
  
  TH1F *hdum = new TH1F("hdum","",1,0.,Rmax*1.1);
  hdum->SetMaximum(1.1*ldf_Ne);
  hdum->SetMinimum(0.1);
  hdum->Draw();
  gStyle->SetOptStat(0);
  TGraphErrors *gLDF = new TGraphErrors(5,ldf_r,ldf_Q,Er,EQ);
  gPad->SetLogy(1);
  gLDF->Draw("*");
  fNKG->Draw("same");

  TLegend *legLDF = new TLegend(0.4,0.7,0.9,0.9);
   legLDF->SetHeader("LDF fit : NKG function, fitting N_{e} and (x,y)_{core}.");
   legLDF->AddEntry((TObject*)0,"NKG(r) = N_{e} #frac{C_{s}}{r_{m}^{2}} #left(#frac{r}{r_{m}}#right)^{s-2} #left(1+#frac{r}{r_{m}}#right)^{s-4.5}","");
   legLDF->AddEntry((TObject*)0,"with s = 1.2","");
   legLDF->AddEntry((TObject*)0,"     C_{s}=0.366 s^{2}(2.07-s)^{1.25} = 0.443","");
   legLDF->AddEntry((TObject*)0,"     r_{m} = 80 m","");
   char sl[80];
   sprintf(sl,"N_{e} = %4.1e #pm %4.1e",ldf_Ne,ldf_ENe);
   legLDF->AddEntry((TObject*)0,sl,"");
   legLDF->Draw();}