//----------------------------------------------------------------------------
// File and Version Information: 
// $Id: TSvtAnaModule.cc,v 1.8 2003/03/21 21:49:49 murat Exp $/
//
// SVTD bank analysis module. Here is what one can do in an interactive session
//
// ---- To run
// TStnAna x("stnmaker.root")
// .L TSvtAnaModule.cc+
// x.AddModule(new TSvtAnaModule)
// x.Run()
//
// ---- To display a histogram
// TSvtAnaModule::Hist_t *c
// TSvtAnaModule *m = (TSvtAnaModule *) x.GetModule("SvtAna")
// c = m->GetHist()
// c->fd0H->Draw()
//
// ---- To process one event
// x.ProcessEntry(1) or x.ProcessEvent(run, evt)  *** not working yet ***
// TSvtAnaModule *m = (TSvtAnaModule *) x.GetModule("SvtAna")
// TSvtDataBlock *svtd = m->GetSvtDataBlock()
// svtd->Print()
//
// Author List:
//    Subir Sarkar  
//----------------------------------------------------------------------------

#include "TCanvas.h"
#include "TF1.h"

#include "Stntuple/obj/TStnTriggerBlock.hh"
#include "Stntuple/obj/TStnHeaderBlock.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"
#include "TSvtAnaModule.hh"
                            // *****  Constructor and destructor
TSvtAnaModule::TSvtAnaModule(const char* name, const char* title):
  TStnModule(name, title)
{
}
TSvtAnaModule::~TSvtAnaModule() {
}
                            // *****  Book Histograms here
void TSvtAnaModule::BookHistograms() {
  Delete("hist");

  // num of tracks 
  fHist.fntrackH = new TH1F("SVTD_ntracks", "SVTD: Number of Tracks ", 10, 0., 10.);
  AddHistogram(fHist.fntrackH);

  // Scatter plot d0 vs phi0
  fHist.fd0phi0H = 
    new TH2F("SVTD_d0phi0", "SVTD: Track d0 vs phi0 ", 50, 0., 2*M_PI, 200, -0.5, 0.5);
  AddHistogram(fHist.fd0phi0H);

  // Scatter plot chi2 vs phi0
  fHist.fchi2phi0H = 
    new TH2F("SVTD_chi2phi0","SVTD: tracks chi2 vs. phi0 ", 50, 0., 2*M_PI, 1000, 0., 1000.); 
  AddHistogram(fHist.fchi2phi0H);

  // Scatter plot d0 vs chi2 
  fHist.fd0chi2H = 
    new TH2F("SVTD_d0chi2","SVTD: tracks d0 vs chi2 ", 1000, 0., 1000., 200, -0.5, 0.5); 
  AddHistogram(fHist.fd0chi2H);

  // Scatter plot d0 vs zin
  fHist.fd0zinH = 
    new TH2F("SVTD_d0zin","SVTD: tracks d0 vs zin ", 6, 0., 6., 200, -0.5, 0.5); 
  AddHistogram(fHist.fd0zinH);

  // Scatter plot d0 vs zout
  fHist.fd0zoutH = 
    new TH2F("SVTD_d0zout","SVTD: tracks d0 vs zout ", 6, 0., 6., 200, -0.5, 0.5); 
  AddHistogram(fHist.fd0zoutH);

  // Wedge number
  fHist.fwedgeH = new TH1F("SVTD_wedges","SVTD: Wedge Occupancy ", 12, 0., 12.);
  AddHistogram(fHist.fwedgeH);

  // Road number 
  fHist.froadH = new  TH1F("SVTD_roads", "SVTD: ROAD Occupancy ", 1000, 0., 40000.);
  AddHistogram(fHist.froadH);

  // Barrel Occupancy (in)
  fHist.fzinH = new TH1F("SVTD_zin","SVTD: Barrel Occupancy ", 6, 0., 6.);
  AddHistogram(fHist.fzinH);

  // Barrel Occupancy (out)
  fHist.fzoutH = new TH1F("SVTD_zout","SVTD: Barrel Occupancy ", 6, 0., 6.);
  AddHistogram(fHist.fzoutH);

  // Impact paramater
  fHist.fd0H = new TH1F("SVTD_d0","SVTD: track Impact Parameter ",100, -0.5, 0.5);
  AddHistogram(fHist.fd0H);

  // phi
  fHist.fphi0H = new TH1F("SVTD_phi0", "SVTD: track phi0 ", 50, 0., 2*M_PI); 
  AddHistogram(fHist.fphi0H);

  // Transverse momentum
  fHist.fptH = new TH1F("SVTD_pt", "SVTD: track Pt ", 100, 0., 10.);
  AddHistogram(fHist.fptH);

  // Curvature 
  fHist.fcurvH = new TH1F("SVTD_curv", "SVTD: track curvature ", 100, -0.003, 0.003);
  AddHistogram(fHist.fcurvH);

  // Fit chi square
  fHist.fchi2H = new TH1F("SVTD_chi2", "SVTD: track Chi2 ", 1000, 0., 1000.);
  AddHistogram(fHist.fchi2H);

  fHist.fxftH = new TH1F("SVTD_xft", "SVTD: XFT wedge number", 300, 0., 300.); // 288 phi slices
  AddHistogram(fHist.fxftH);

  fHist.fhitH[0] = new TH1F("SVTD_hit0", "SVTD: Hit 0 ", 100, 0., 0.1);
  AddHistogram(fHist.fhitH[0]);

  fHist.fhitH[1] = new TH1F("SVTD_hit1", "SVTD: Hit 1 ", 100, 0., 0.1);
  AddHistogram(fHist.fhitH[1]);

  fHist.fhitH[2] = new TH1F("SVTD_hit2", "SVTD: Hit 2 ", 100, 0., 0.1);
  AddHistogram(fHist.fhitH[2]);

  fHist.fhitH[3] = new TH1F("SVTD_hit3", "SVTD: Hit 3 ", 100, 0., 0.1);
  AddHistogram(fHist.fhitH[3]);

  fHist.fhitH[4] = new TH1F("SVTD_hit4", "SVTD: Hit 4 ", 100, 0., 0.1);
  AddHistogram(fHist.fhitH[4]);

  // Hit positions with respect to wedge border
  fHist.fwhitH[0] = new TH1F("SVTD_whit0", "SVTD: Hit 0 ", 200, 0., 2.);
  AddHistogram(fHist.fwhitH[0]);

  fHist.fwhitH[1] = new TH1F("SVTD_whit1", "SVTD: Hit 1 ", 250, 0., 2.5);
  AddHistogram(fHist.fwhitH[1]);

  fHist.fwhitH[2] = new TH1F("SVTD_whit2", "SVTD: Hit 2 ", 400, 0., 4.0);
  AddHistogram(fHist.fwhitH[2]);

  fHist.fwhitH[3] = new TH1F("SVTD_whit3", "SVTD: Hit 3 ", 500, 0., 5.0); 
  AddHistogram(fHist.fwhitH[3]);

  fHist.fwhitH[4] = new TH1F("SVTD_whit4", "SVTD: Hit 4 ", 600, 0., 6.0); 
  AddHistogram(fHist.fwhitH[4]);

  // Fit Quality from TF status word
  fHist.fFitQH = new TH1F("SVTD_FitQuality", "SVTD: Fit Quality", 16, 0., 16.);
  AddHistogram(fHist.fFitQH);

  // TF error word
  fHist.fTFerrH = new TH1F("SVTD_TFerr", "SVTD: TF error word", 7, 0., 7.); 
  AddHistogram(fHist.fTFerrH);

  // End Event error word
  fHist.fEEerrH = new TH1F("SVTD_EEerr", "SVTD: End Event error word", 8, 0., 8.);
  AddHistogram(fHist.fEEerrH);

                                 // ****** XFT quantities
  // Number of tracks
  fHist.fntXftH = new TH1F("SVTD_ntracks_xft", "SVTD: Number of XFT tracks", 50, 0., 50.);
  AddHistogram(fHist.fntXftH);

  // Phi bin
  fHist.fphibinH = new TH1F("SVTD_phibin_xft", "SVTD: Phi bin for XFT track", 100, 0., 4500.);
  AddHistogram(fHist.fphibinH);

  // Mini Wedge
  fHist.fminiwedgeH = new TH1F("SVTD_mwedge_xft", "SVTD: Mini wegde for XFT track", 100, 0., 600.);
  AddHistogram(fHist.fminiwedgeH);

  // Mini phi
  fHist.fminiphiH = new TH1F("SVTD_mphi_xft", "SVTD: Mini phi for XFT track", 8, 0., 8.);
  AddHistogram(fHist.fminiphiH);

  // Wedge
  fHist.fwedgeXftH = new TH1F("SVTD_wedge_xft", "SVTD: wegde for XFT track", 12, 0., 12.);
  AddHistogram(fHist.fwedgeXftH);

  // Pt bin
  fHist.fptbinH = new TH1F("SVTD_ptbin_xft", "SVTD: Pt bin for XFT track", 100, 0., 100.);
  AddHistogram(fHist.fptbinH);

  // Isolation flag
  fHist.fisoH = new TH1F("SVTD_iso_xft", "SVTD: Isolation flag for XFT track", 2, 0., 2.);
  AddHistogram(fHist.fisoH);

  // Short track flag
  fHist.fshtrkH = new TH1F("SVTD_shtrk_xft", "SVTD: Short track flag for XFT track", 2, 0., 2.);
  AddHistogram(fHist.fshtrkH);

  // Phi in radian
  fHist.fphiradH = new TH1F("SVTD_phirad_xft", "SVTD: Phi for XFT track", 100, 0., 2.*M_PI);
  AddHistogram(fHist.fphiradH);

  // Curvature
  fHist.fcurvXftH = new TH1F("SVTD_curv_xft", "SVTD: Curv for XFT track", 100, -0.0015, 0.0015);
  AddHistogram(fHist.fcurvXftH);

  // Track pt
  fHist.fptXftH = new TH1F("SVTD_pt_xft", "SVTD: Pt for XFT track", 100, -20., 20.);
  AddHistogram(fHist.fptXftH);
}

//_____________________________________________________________________________
int TSvtAnaModule::BeginJob() {
  RegisterDataBlock("SvtDataBlock","TSvtDataBlock",&fSvtData);
                                // book histograms
  BookHistograms();
  if (!fSvtData) {
    printf(" --> Branch *** TSvtDataBlock *** doesn't exist!\n");
    fEnabled = 0;
  }
  return 0;
}

//_____________________________________________________________________________
int TSvtAnaModule::BeginRun() {
  return 0;
}

//_____________________________________________________________________________
int TSvtAnaModule::Event(int ientry) {
  fSvtData->GetEntry(ientry);
  if (PrintLevel()) {
    GetHeaderBlock()->Print();
    fSvtData->Print();
  }
                                       // ***** SVT tracks
  int nSvtTracks = fSvtData->NSvtTracks();
  fHist.fntrackH->Fill(nSvtTracks*1.0);    
  for (int i = 0; i < nSvtTracks; i++) 
  {
    TSvtTrack *trk = fSvtData->SvtTrack(i);
    Float_t d0   = trk->D0();
    Float_t phi0 = trk->Phi();
    Float_t chi2 = trk->Chi2();
    Float_t curv = trk->Curv();
    Float_t pt   = trk->Pt();
    Int_t wedge  = trk->Wedge();
    Int_t road   = trk->Road();
    Int_t zin    = trk->Zin();
    Int_t zout   = trk->Zout();
    Int_t xft    = trk->XftNumber();
    Int_t hits[5];
    Float_t hitx[5];
    for (int j = 0; j < 5; j++) {
      hits[j] = trk->Hit16(j);
      fHist.fwhitH[j]->Fill(hits[j]);     // Hit position j wrt wedge border!

      hitx[j] = trk->Hit(j);  
      fHist.fhitH[j]->Fill(hitx[j]);      // Hit Position j
    }
    Int_t fit_q  = trk->FitQuality();
    Int_t TF_err = trk->TFError();

    fHist.fd0H->Fill(d0+0.00001);         // Distr. of impact parameters
    fHist.fphi0H->Fill(phi0);             // Phi
    fHist.fd0phi0H->Fill(phi0, d0);       // Scatter plot d0 vs phi0 
    fHist.fchi2phi0H->Fill(phi0, chi2);   // Scatter plot chi2 vs phi0
    fHist.fd0chi2H->Fill(chi2, d0);       // Scatter plot d0 vs chi2
    fHist.fd0zinH->Fill(zin*1.0, d0);     // Scatter plot d0 vs zin
    fHist.fd0zoutH->Fill(zout*1.0, d0);   // Scatter plot d0 vs zout
    fHist.fchi2H->Fill(chi2);             // Fit Chi Square    
    fHist.fwedgeH->Fill(wedge*1.0);       // occupancy per wedge
    fHist.froadH->Fill(road*1.0);         // Road ID
    fHist.fzinH->Fill(zin*1.0);           // Barrel occupancy (in)
    fHist.fzoutH->Fill(zout*1.0);         // Barrel occupancy (out)
    fHist.fptH->Fill(pt);                 // Pt
    fHist.fcurvH->Fill(curv);             // Curvature
    fHist.fxftH->Fill(xft*1.0);           // XFT track number
    fHist.fFitQH->Fill(fit_q*1.0);        // Fit quality from TF status word

    // Error bits in TF status word
    for (int j = 0; j < 7; j++) {
      if ((TF_err >> j) & 0x1) fHist.fTFerrH->Fill(j*1.0);
    }
  }

  // End Event Error
  Int_t EEerr = fSvtData->GError();
  for (int i = 0; i < 8; i++) {
    if ((EEerr >> i) & 0x1) fHist.fEEerrH->Fill(i*1.0);
  }
                                      // ***** XFT tracks
  int nXftTracks = fSvtData->NXftTracks();
  int nTracks = 0;
  for (int i = 0; i < nXftTracks; i++) 
  {
    TSvXtTrack *trk = fSvtData->XftTrack(i);
    if (!trk->IsValid()) continue;
    nTracks++;
    Int_t phi       = trk->Phi();
    Int_t miniwedge = trk->Miniwedge();
    Int_t miniphi   = trk->Miniphi();
    Int_t wedge     = trk->Wedge();
    Int_t ptbin     = trk->Ptbin();
    Int_t iso       = (trk->IsIsolated()) ? 1 : 0;
    Int_t shtrk     = (trk->IsShort()) ? 1 : 0;
    Float_t phirad  = trk->Phirad();
    Float_t curv    = trk->Curv();
    Float_t pt      = trk->Pt();

    fHist.fphibinH->Fill(phi*1.0);           // Phi bin
    fHist.fminiwedgeH->Fill(miniwedge*1.0);  // Mini wedge
    fHist.fminiphiH->Fill(miniphi*1.0);      // Mini phi
    fHist.fwedgeXftH->Fill(wedge*1.0);       // Wedge Number
    fHist.fptbinH->Fill(ptbin*1.0);          // Pt bin
    fHist.fisoH->Fill(iso*1.0);              // isolation flag
    fHist.fshtrkH->Fill(shtrk*1.0);          // Short track flag
    fHist.fphiradH->Fill(phirad);            // phi in radian
    fHist.fcurvXftH->Fill(curv);             // Curvature
    fHist.fptXftH->Fill(pt);                 // Transverse Momentum
  }
  fHist.fntXftH->Fill(nTracks*1.0);    

  return 0;
}
                                   // ****  Display histograms
void TSvtAnaModule::DisplayEvent() {
}
                                   // ****  Bookkeeping
int TSvtAnaModule::EndJob() {
  printf("----- End job: ---- %s\n",GetName());
  return 0;
}