//_____________________________________________________________________________
// track-based conversion filter module
// This is an extraction from Ray's TPhotonAnaModule placed as a separate 
// module to facilitate optimization of the cuts. 
// May be we can merge it back?
// The module selects conversions as pairs of opposite-sign tracks which
// are quasi-parallel in the intersection point. In more detail:
// - each track is required to have 2 or more axial and 2 or more stereo
//   COT segments (whatever the definition of the segment is)
// - |ctg(theta1)-ctg(theta2)| < fMaxDeltaLam (default (0.1)
// - dxy < fMaxSeparation, where dxy is distance between the tracks in the
//   point where they are parallel in XY (default 5mm)
//_____________________________________________________________________________

#include "TF1.h"
#include "TCanvas.h"
#include "Stntuple/loop/TStnAna.hh"
#include "Stntuple/obj/TStnTrackBlock.hh"

#include "TConversionFilterModule.hh"

// ClassImp(TConversionFilterModule)
//_____________________________________________________________________________
TConversionFilterModule::TConversionFilterModule(const char* name, const char* title):
  TStnModule(name,title)
{
  fMaxSeparation = 0.5;
  fMaxDeltaLam   = 0.1;
}

//_____________________________________________________________________________
TConversionFilterModule::~TConversionFilterModule() {
}

//_____________________________________________________________________________
void TConversionFilterModule::BookHistograms() {
  char name [200];
  char title[200];

  Delete("hist");

  sprintf(name, "%s_nst_vs_nax",GetName());
  sprintf(title,"%s: COT track: N(st) vs N(ax) segments",GetName());
  fHist.fNSeg = new TH2F(name, title,5,0,5,5,0,5);
  AddHistogram(fHist.fNSeg);

  sprintf(name, "%s_dlam_rs_0",GetName());
  sprintf(title,"%s: delta(lambda), charge = 0, (0)",GetName());
  fHist.fDeltaLamRS[0] = new TH1F(name, title,250,-2.5,2.5);
  AddHistogram(fHist.fDeltaLamRS[0]);

  sprintf(name, "%s_dlam_ws_0",GetName());
  sprintf(title,"%s: delta(lambda), charge != 0, (0)",GetName());
  fHist.fDeltaLamWS[0] = new TH1F(name, title,250,-2.5,2.5);
  AddHistogram(fHist.fDeltaLamWS[0]);

  sprintf(name, "%s_dlam_rs_1",GetName());
  sprintf(title,"%s: delta(lambda), charge = 0, (1)",GetName());
  fHist.fDeltaLamRS[1] = new TH1F(name, title,100,-0.25,0.25);
  AddHistogram(fHist.fDeltaLamRS[1]);

  sprintf(name, "%s_dlam_ws_1",GetName());
  sprintf(title,"%s: delta(lambda), charge != 0, (1)",GetName());
  fHist.fDeltaLamWS[1] = new TH1F(name, title,100,-0.25,0.25);
  AddHistogram(fHist.fDeltaLamWS[1]);

  sprintf(name, "%s_sep_rs_0",GetName());
  sprintf(title,"%s: conv cand: separation /right sign/0",GetName());
  fHist.fSeparationRS[0] = new TH1F(name, title,250,-25,25);
  AddHistogram(fHist.fSeparationRS[0]);

  sprintf(name, "%s_sep_ws_0",GetName());
  sprintf(title,"%s: conv cand: separation /wrong sign/0",GetName());
  fHist.fSeparationWS[0] = new TH1F(name, title,250,-25,25);
  AddHistogram(fHist.fSeparationWS[0]);

  sprintf(name, "%s_sep_rs_1",GetName());
  sprintf(title,"%s: conv cand: separation /right sign/1",GetName());
  fHist.fSeparationRS[1] = new TH1F(name, title,100,-5,5);
  AddHistogram(fHist.fSeparationRS[1]);

  sprintf(name, "%s_sep_ws_1",GetName());
  sprintf(title,"%s: conv cand: separation /wrong sign/0",GetName());
  fHist.fSeparationWS[1] = new TH1F(name, title,100,-5,5);
  AddHistogram(fHist.fSeparationWS[1]);

  for (int i=0; i<7; i++) {
    sprintf(name, "%s_rconv_rs_%i",GetName(),i);
    sprintf(title,"%s: Conversion radius (%i) right sign",GetName(),i);
    fHist.fRConvRS[i] = new TH1F(name, title,500,0,100);
    AddHistogram(fHist.fRConvRS[i]);

    sprintf(name, "%s_rconv_ws_%i",GetName(),i);
    sprintf(title,"%s: Conversion radius (%i) wrong sign",GetName(),i);
    fHist.fRConvWS[i] = new TH1F(name, title,500,0,100);
    AddHistogram(fHist.fRConvWS[i]);

    sprintf(name, "%s_ncand_%i",GetName(),i);
    sprintf(title,"%s: N(conv candidates) (%i)",GetName(),i);
    fHist.fNCand[i] = new TH1F(name, title,100,0,100);
    AddHistogram(fHist.fNCand[i]);

  }

  sprintf(name, "%s_yconv_vs_xconv",GetName());
  sprintf(title,"%s: Conversion Y:X",GetName());
  fHist.fYConvVsXConv = new TH2F(name, title,100,-50,50,100,-50,50);
  AddHistogram(fHist.fYConvVsXConv);


}


//_____________________________________________________________________________
int TConversionFilterModule::BeginJob() {

				// register the data block

  fTrackBlock   = (TStnTrackBlock*) RegisterDataBlock("TrackBlock",
						      "TStnTrackBlock");

  fTriggerBlock = (TStnTriggerBlock*) RegisterDataBlock("TriggerBlock",
							"TStnTriggerBlock");

					// book histograms
  BookHistograms();

  return 0;
}


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

//_____________________________________________________________________________
int TConversionFilterModule::Event(int ientry) {
  // if fMyronFlag >= 0, select event with the given L2 bucket = fMyronFlag

  if (fMyronFlag >= 0) {
    fTriggerBlock->GetEntry(ientry);
    if (fTriggerBlock->Tsid()->MyronFlag() != fMyronFlag) {
      return 0;
    }
  }
    

  fTrackBlock->GetEntry(ientry);

  double  r1, r2, x1, x2, y1, y2, ff, dd, sep, dt, conv_sign;
  double  xconv, yconv, rconv, dx, dy, sqr;
  int     ncand, ntrk;

  ncand = 0;

  ntrk = fTrackBlock->NTracks();

  for(int i=0; i<ntrk; i++) {
    TStnTrack* t = fTrackBlock->Track(i);
    fHist.fNSeg->Fill(t->NStSeg(),t->NAxSeg(),1);
  }

  for(int i=0; i<ntrk-1; i++) {
    TStnTrack& t1 = *fTrackBlock->Track(i);

    if ((t1.NAxSeg() <= 1) || (t1.NStSeg() <= 1)) goto NEXT_T1;

					// loop over the 2nd track
    for(int j=i+1; j<ntrk; j++) {
      TStnTrack& t2 = *fTrackBlock->Track(j);

      if ((t2.NAxSeg() <= 1) || (t2.NStSeg() <= 1)) goto NEXT_T2;

      dt        = t1.Lam0()-t2.Lam0();
      conv_sign = t1.Charge()+t2.Charge();

      if (conv_sign == 0) fHist.fDeltaLamRS[0]->Fill(dt);
      else                fHist.fDeltaLamWS[0]->Fill(dt);

      if (fabs(dt) < 0.25) {

	r1 = 1.0/fabs(2*t1.C0());
	ff = t1.Charge()*r1+t1.D0();
	x1 = -ff*sin(t1.Phi0());
	y1 =  ff*cos(t1.Phi0());

	r2 = 1.0/fabs(2*t2.C0());
	ff = t2.Charge()*r2+t2.D0();
	x2 = -ff*sin(t2.Phi0());
	y2 =  ff*cos(t2.Phi0());

	dd  = sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));

	sep = dd-r1-r2;

	dx  = x2-x1;
	dy  = y2-y1;
	sqr = sqrt(dx*dx+dy*dy);

	xconv = x1+dx/sqr*(r1+sep/2);
	yconv = y1+dy/sqr*(r1+sep/2);

	rconv = sqrt(xconv*xconv+yconv*yconv);

	if (conv_sign == 0) fHist.fSeparationRS[0]->Fill(sep);
	else                fHist.fSeparationWS[0]->Fill(sep);
//-----------------------------------------------------------------------------
// take some "arbitrary" cut on delta(lambda) and plot distribution for 
// separation in XY
//-----------------------------------------------------------------------------
	if (fabs(dt) < 0.1) {
	  if (conv_sign == 0) fHist.fSeparationRS[1]->Fill(sep);
	  else                fHist.fSeparationWS[1]->Fill(sep);
	}
//-----------------------------------------------------------------------------
// now for the track pairs with small dR(xy) plot distribution for
// delta(Lambda)
//-----------------------------------------------------------------------------
	if (fabs(sep) < 0.5) {
	  if (conv_sign == 0) fHist.fDeltaLamRS[1]->Fill(dt);
	  else                fHist.fDeltaLamWS[1]->Fill(dt);
	}
//-----------------------------------------------------------------------------
// count number of candidate conversion pairs per event
//-----------------------------------------------------------------------------
	if ( (fabs(sep) < fMaxSeparation) && (fabs(dt) < fMaxDeltaLam)) {
	  ncand += 1;
	  if (conv_sign == 0) fHist.fRConvRS[0]->Fill(rconv);
	  else                fHist.fRConvWS[0]->Fill(rconv);
	}
      }
    NEXT_T2:;
    }
  NEXT_T1:;
  }

  for (int i=0; i<7; i++) {
    fHist.fNCand[i]->Fill(ncand);
  }

  if (ncand > 0) SetPassed(1);
  else           SetPassed(0);

  return 0;
}


//_____________________________________________________________________________
int TConversionFilterModule::EndJob() {
  return 0;
}