////////////////////////////////////////////////////////////////
//
//  Histogram set for tracker-efficiency study.
//
//  It contains the functions:
//
//  - CreateHistos()
//  - FillHistos(PamLevel2* event)
//  - Save Histos()
//
//  Usage (example):
//  [] .L My-histo-set-macro.C+
//  [] .L My-selection.C+
//  [] .L Loop.C+
//  [] Loop("data-dir/","file-list.txt",100000,"+ALL fillTree fillHistos","output-file.root")
//  
////////////////////////////////////////////////////////////////

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMath.h>
#include <TLine.h>
#include <TPolyMarker.h>
#include <TSelector.h>
#include <TFile.h>

#include <stdlib.h>
#include <iostream>
using namespace std;

#include <PamLevel2.h>
#include <CaloAxis.h>

#endif
//======================
// HITOGRAMS declaration
//======================
TH1F *htofresx;
TH1F *htofresy;
TH1F *hcalresx;
TH1F *hcalresy;
TH1F *htrkresx[10];
TH1F *htrkresy[10];
TH1F *hxcalotrk[6];
TH1F *hycalotrk[6];
TH1F *hresangx;
TH1F *hresangy;
//===============================================================
// Create histograms
//===============================================================
void CreateHistos(  ){

    gROOT->cd();//create histos in memory

    gDirectory->Delete("htofresx");
    gDirectory->Delete("htofresy");
    TH1F *htofresx = new TH1F("htofresx","Spatial residuals (S12+S32) - S21 (x)",200,-10.,10.);
    TH1F *htofresy = new TH1F("htofresy","Spatial residuals (S11+S31) - S22 (y)",200,-10.,10.);


    gDirectory->Delete("hcalresx");
    gDirectory->Delete("hcalresy");
    TH1F *hcalresx = new TH1F("hcalresx","Spatial residuals on CAL - S32 (x)",200,-10.,10.);
    TH1F *hcalresy = new TH1F("hcalresy","Spatial residuals on CAL - S31 (y)",200,-10.,10.);

    TString hid; 
    for(Int_t i=0; i<3; i++){
	hid.Form("htrkresx%i",i);
	htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - TOF (x)",200,-10.,10.);
	hid.Form("htrkresy%i",i);
	htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - TOF (y)",200,-10.,10.);	
    }
    for(Int_t i=3; i<9; i++){
	hid.Form("htrkresx%i",i);
	htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK (x)",5000,-.2,.2);
	hid.Form("htrkresy%i",i);
	htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK (y)",5000,-.2,.2);	
    }
    for(Int_t i=9; i<10; i++){
	hid.Form("htrkresx%i",i);
	htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - CAL (x)",800,-10.,10.);
	hid.Form("htrkresy%i",i);
	htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - CAL (y)",800,-10.,10.);	
    }


    for(Int_t ip=0; ip<6; ip++){
	hid.Form("hxcalotrk%i",ip);
	hxcalotrk[ip] = new TH1F(hid.Data(),"Calo axis intersaction on Trk plane (x)",500,-25.,25.);
	hid.Form("hycalotrk%i",ip);
	hycalotrk[ip] = new TH1F(hid.Data(),"Calo axis intersaction on Trk plane (y)",500,-25.,25.);
    }

    hresangx = new TH1F("hresangx","Angular residual on the top calo plane (x)",400,-40.,40.);
    hresangy = new TH1F("hresangy","Angular residual on the top calo plane (y)",400,-40.,40.);

}

//===============================================================
// Fill histograms (called event-by-event)
//===============================================================
 
void FillHistos( PamLevel2* event ){
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    htofresx = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresx"));
    htofresy = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresy"));
//    cout << hex << htofresx << dec << endl;
    hcalresx = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresx"));
    hcalresy = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresy"));
    TString hid; 
    for(Int_t i=0; i<10; i++){
	hid.Form("htrkresx%i",i);
	htrkresx[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
	hid.Form("htrkresy%i",i);
	htrkresy[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));	
    }
    for(Int_t ip=0; ip<6; ip++){
	hid.Form("hxcalotrk%i",ip);
	hxcalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
	hid.Form("hycalotrk%i",ip);
	hycalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
    }
    hresangx = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangx"));
    hresangy = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangy"));

//  ----------------------------------
//  calo variables
//  ----------------------------------
    CaloAxis *x_axis = new CaloAxis();
    CaloAxis *y_axis = new CaloAxis();

    float rcil = 1.;// tolerance (cm)
    x_axis->FitAxis(event->GetCaloLevel1(),0,rcil);
    y_axis->FitAxis(event->GetCaloLevel1(),1,rcil);

    float xs31 = x_axis->par[0]+x_axis->par[1]*event->GetZTOF(31);
    float xs32 = x_axis->par[0]+x_axis->par[1]*event->GetZTOF(32);
    float ys31 = y_axis->par[0]+y_axis->par[1]*event->GetZTOF(31);
    float ys32 = y_axis->par[0]+y_axis->par[1]*event->GetZTOF(32);
    float qtrack = x_axis->GetQaxis()+y_axis->GetQaxis();
    int   ntrack = x_axis->GetN()+y_axis->GetN();

    CaloStrip top = CaloStrip();
    top.Set(0,0,48);
    float zcalotop = top.GetZ();

    float xcalotop = x_axis->par[0]+x_axis->par[1]*zcalotop;
    float ycalotop = y_axis->par[0]+y_axis->par[1]*zcalotop;
    
    
    float xtrk[6];
    float ytrk[6];
    for(Int_t ip=0; ip<6; ip++){
	xtrk[ip] =  x_axis->par[0]+x_axis->par[1]*event->GetZTrk(ip+1);
	ytrk[ip] =  y_axis->par[0]+y_axis->par[1]*event->GetZTrk(ip+1);
    }

    float caloangx = (float)atan((double)x_axis->par[1])*180./3.1415026;
    float caloangy = (float)atan((double)y_axis->par[1])*180./3.1415026;

    float xtrktop = x_axis->par[0]+x_axis->par[1]*event->GetZTrk(1);
    float ytrktop = y_axis->par[0]+y_axis->par[1]*event->GetZTrk(1);

    x_axis->Delete();
    y_axis->Delete();

    // tolleranza sul contenimento della traccia nel tracker
    float dztrk = 44.51;
    float dxtrk = 16.14;
    float dytrk = 13.14;
    float axtoll = 10.; //deg
    float aytoll = 5.;  //deg
    float tollx = tan(atan(dxtrk/dztrk)+axtoll*3.1415026/180.)*dztrk-dxtrk;
    float tolly = tan(atan(dytrk/dztrk)+aytoll*3.1415026/180.)*dztrk-dytrk;
//    cout << tollx << " "<<tolly << endl;

//  ----------------------------------
//  tof variables
//  ----------------------------------
//  get track stored by the tof    
    ToFTrkVar *tof = event->GetToFStoredTrack(-1);
    if(!tof){
	cout << " no ToF-track stored "<<endl;
	return;
    }
    float resxs21 = (event->GetZTOF(21)-event->GetZTOF(32))*(tof->xtofpos[0]-tof->xtofpos[2])/(event->GetZTOF(12)-event->GetZTOF(32))+tof->xtofpos[2];
    float resys22 = (event->GetZTOF(22)-event->GetZTOF(31))*(tof->ytofpos[0]-tof->ytofpos[2])/(event->GetZTOF(11)-event->GetZTOF(31))+tof->ytofpos[2];
    resxs21 -= tof->xtofpos[1];
    resys22 -= tof->ytofpos[1];

    float resxs32 = xs32 - tof->xtofpos[2];
    float resys31 = ys31 - tof->ytofpos[2];

//  ----------------------------------
//  tracker variables
//  ----------------------------------
    float resxtrk[]={-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.};
    float resytrk[]={-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.};

    float resangx = -999;
    float resangy = -999;
    if( event->GetNTracks()==1 ){
	
	PamTrack *track = event->GetTrack(0);
	if( track->chi2 >0  ){

	    Float_t ztraj[13];
	    
	    Int_t i=0;
	    ztraj[i++] = event->GetZTOF(11);
	    ztraj[i++] = event->GetZTOF(12);
	    ztraj[i++] = event->GetZTOF(21);
	    ztraj[i++] = event->GetZTOF(22);
	    ztraj[i++] = event->GetZTrk(1);
	    ztraj[i++] = event->GetZTrk(2);
	    ztraj[i++] = event->GetZTrk(3);
	    ztraj[i++] = event->GetZTrk(4);
	    ztraj[i++] = event->GetZTrk(5);
	    ztraj[i++] = event->GetZTrk(6);
	    ztraj[i++] = event->GetZTOF(31);
	    ztraj[i++] = event->GetZTOF(32);
	    ztraj[i++] = zcalotop;
	    Trajectory traj = Trajectory(13,ztraj);
	    traj.DoTrack2( track->al );
	    
	    // --- S1
	    resxtrk[0] = track->xtofpos[0] - traj.x[1];
	    resytrk[0] = track->ytofpos[0] - traj.y[0];
	    // --- S2
	    resxtrk[1] = track->xtofpos[1] - traj.x[2];
	    resytrk[1] = track->ytofpos[1] - traj.y[3];
	    // --- S3
	    resxtrk[2] = track->xtofpos[2] - traj.x[11];
	    resytrk[2] = track->ytofpos[2] - traj.y[10];
	    // --- trk
	    for(Int_t ip=0; ip<6; ip++){
		if(track->XGood(ip))resxtrk[3+ip] = track->xm[ip] - traj.x[4+ip];
		if(track->YGood(ip))resytrk[3+ip] = track->ym[ip] - traj.y[4+ip];
	    }
	    // --- calo
	    resxtrk[9] = xcalotop - traj.x[12];
	    resytrk[9] = ycalotop - traj.y[12];

	    resangx = caloangx - traj.thx[12];
	    resangy = caloangy - traj.thy[12];
//	    cout << caloangx << " -x- " << traj.thx[12]<<endl;
//	    cout << caloangy << " -y- " << traj.thy[12]<<endl;

	}
	if(track) delete track;
    }
    
//  ----------------------------------
//  fill spatial residuals 
//  ----------------------------------
    htofresx->Fill(resxs21);
    htofresy->Fill(resys22);
    hcalresx->Fill(resxs32);
    hcalresy->Fill(resys31);
    for(Int_t ih=0; ih<10; ih++){
	htrkresx[ih]->Fill(resxtrk[ih]);
	htrkresy[ih]->Fill(resytrk[ih]);
    }
    for(Int_t ip=0; ip<6; ip++){
	hxcalotrk[ip]->Fill(xtrk[ip]);
	hycalotrk[ip]->Fill(ytrk[ip]);
    }
    hresangx->Fill(resangx);
    hresangy->Fill(resangy);


}
//===============================================================
// Save histograms
//===============================================================
void SaveHistos( TFile * outf ){

    TString hid; 

    gROOT->cd();
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    htofresx = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresx"));
    htofresy = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresy"));
//    cout << hex << htofresx << dec << endl;
    hcalresx = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresx"));
    hcalresy = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresy"));
    for(Int_t i=0; i<10; i++){
	hid.Form("htrkresx%i",i);
	htrkresx[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
	hid.Form("htrkresy%i",i);
	htrkresy[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));	
    }
    for(Int_t ip=0; ip<6; ip++){
	hid.Form("hxcalotrk%i",ip);
	hxcalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
	hid.Form("hycalotrk%i",ip);
	hycalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
    }
    hresangx = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangx"));
    hresangy = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangy"));

//  ----------------------------------
//  save on file
//  ----------------------------------

    outf->cd();
    
    if(htofresx)htofresx->Write();
    if(htofresy)htofresy->Write();
    if(hcalresx)hcalresx->Write();
    if(hcalresy)hcalresy->Write();

    for(Int_t i=0; i<10; i++){
	if(htrkresx[i])htrkresx[i]->Write();
	if(htrkresy[i])htrkresy[i]->Write();
    }

    for(Int_t ip=0; ip<6; ip++){
	if(hxcalotrk[ip])hxcalotrk[ip]->Write();
	if(hycalotrk[ip])hycalotrk[ip]->Write();
    }

    if(hresangx) hresangx->Write();
    if(hresangy) hresangy->Write();
}