//_____________________________________________________________________________
// SVX analysis module
//_____________________________________________________________________________
#include "TF1.h"
#include "Stntuple/loop/TStnAna.hh"

#include "TSvxAnaModule.hh"
#include "Stntuple/obj/TStnHeaderBlock.hh"

#include <iostream>
#include <cmath>
#include <algorithm>

//_____________________________________________________________________________
TSvxAnaModule::TSvxAnaModule(const char* name, const char* title):
  TStnModule(name,title), fFile(NULL), fAlg(11), fAreSiHits(1), 
  fTreeA(true), fTreeT(true), fTreeN(true), fTreeR(true)
{
  fHist.fT=0;
  fHist.fA=0;
  fHist.fN=0;
  fHist.fR=0;
}

//_____________________________________________________________________________
TSvxAnaModule::~TSvxAnaModule() {
  if (fHist.fT) delete fHist.fT;
  if (fHist.fA) delete fHist.fA;
  if (fHist.fN) delete fHist.fN;
  if (fHist.fR) delete fHist.fR;
  if (fFile) fFile->Close();
}




//_____________________________________________________________________________
void TSvxAnaModule::BookHistograms() {
  
  if (fTreeT || fTreeA || fTreeN || fTreeR) {
    fFile = new TFile("fred2.root","recreate");
  }

  if (fHist.fEff_trk) delete fHist.fEff_trk;
  fHist.fEff_trk=new TH1F("eff","eff match",8,-0.5,7.5);

  if (fTreeT) {
    if (fHist.fT) delete fHist.fT;
    fHist.fT = new TTree("T","My hits");
    fHist.fT->Branch("hit",&fHit,"run/I:evt:b:h:p:l:s:stp/F:N/I:Q/F:f:x:y:z:phi0:z0:d0:eta:pt");
  }

  if (fTreeA) {
    if (fHist.fA) delete fHist.fA;
    fHist.fA = new TTree("A","All hits");
    fHist.fA->Branch("all",&fAll,"run/I:evt:b:h:p:l:s:stp/F:N/I:Q/F:f:x:y:z:phi0:z0:d0:eta:pt");
  }

  if (fTreeN) {
    if (fHist.fN) delete fHist.fN;
    fHist.fN = new TTree("N","Near hits");
    fHist.fN->Branch("near",&fNear,"run/I:evt:b:h:p:l:s:act:ix/F:iy:iz:istp:ily:ilz:f:phi0:z0:d0:eta:pt:stp[3]:N[3]/I:Q[3]/F:x[3]:y[3]:z[3]");
  }

  if (fTreeR) {
    if (fHist.fR) delete fHist.fR;
    fHist.fR = new TTree("R","Tracks");
    fHist.fR->Branch("track",&fTrack,"run/I:evt:b:C/F:phi0:d0:tanl:z0:pt");
  }
}

void TSvxAnaModule::SetTrackAlgorithm(int alg) {
  fAlg=alg;
}
void TSvxAnaModule::SetAreSiHits(int areSiHits) {
  fAreSiHits=areSiHits;
}
void TSvxAnaModule::SetTreeA(bool makeTreeA) {
  fTreeA=makeTreeA;
}
void TSvxAnaModule::SetTreeT(bool makeTreeT) {
  fTreeT=makeTreeT;
}
void TSvxAnaModule::SetTreeN(bool makeTreeN) {
  fTreeN=makeTreeN;
}
void TSvxAnaModule::SetTreeR(bool makeTreeR) {
  fTreeR=makeTreeR;
}


//_____________________________________________________________________________
int TSvxAnaModule::BeginJob() {

  // register the data blocks
  fSvxDataBlock=(TSvxDataBlock*) fAna->RegisterDataBlock("SvxDataBlock","TSvxDataBlock");
  fTrackBlock=(TStnTrackBlock*)  fAna->RegisterDataBlock("TrackBlock","TStnTrackBlock");
  fTrackLinkBlock=(TStnTrackLinkBlock*)  fAna->RegisterDataBlock("TrackLinkBlock","TStnTrackLinkBlock");
  fSiIsectBlock=(TSiIsectBlock*) fAna->RegisterDataBlock("SiIsectBlock","TSiIsectBlock");

  // book histograms
  BookHistograms();

  if (! fSvxDataBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","SvxDataBlock");
    fEnabled = 0;
  }
  if (! fTrackBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","TrackBlock");
    fEnabled = 0;
  }
  if (! fTrackLinkBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","TrackLinkBlock");
    fEnabled = 0;
  }
  if (! fSiIsectBlock) {
    printf(" >>> branch *** %s *** doesn't exist \n","SiIsectBlock");
    fEnabled = 0;
  }
  return 0;
}


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

//_____________________________________________________________________________
int TSvxAnaModule::Event(int ientry) {

  if (PrintLevel()>0) {
    printf(" ---- Starting ientry %d\n",ientry);
  }

  fSvxDataBlock->GetEntry(ientry);
  int nSiHits = fSvxDataBlock->NSiHits();
  if (nSiHits<=0) return 0;
  // First make tree of all hits
  if (fTreeA) {
    for (int i=0; i<nSiHits; ++i){
      TStnSiHit* hit = fSvxDataBlock->SiHit(i);
      fAll.run=GetHeaderBlock()->RunNumber();
      fAll.evt=GetHeaderBlock()->EventNumber();
      fAll.b=hit->DigiCode()->Barrel();
      fAll.h=hit->DigiCode()->LadderSeg();
      fAll.p=hit->DigiCode()->PhiWedge();
      fAll.l=hit->DigiCode()->Layer();
      fAll.s=hit->DigiCode()->Side();
      fAll.stp=hit->StripNum();
      fAll.N=hit->Nstrip();
      fAll.Q=hit->Qtot()*1.0;
      fAll.f=1.0;
      fAll.x=hit->Global()->X();
      fAll.y=hit->Global()->Y();
      fAll.z=hit->Global()->Z();
      fHist.fA->Fill();
    }    
  }
    
  // Then check on other branches for data
  fTrackBlock->GetEntry(ientry);
  if (fTrackBlock->NTracks()<1) return 0;
  if (PrintLevel()>2) printf("NTracks=%d\n",fTrackBlock->NTracks());

  fTrackLinkBlock->GetEntry(ientry);
  if (fAreSiHits) {
    if (fTrackLinkBlock->SiHitLinkTrk()->NLinksTotal()<1) return 0;
  }
  if (fTrackLinkBlock->SiIsectLinkTrk()->NLinksTotal()<1) return 0;
  if (PrintLevel()>2) 
    printf("NHitLinks=%d\n",fTrackLinkBlock->SiHitLinkTrk()->NLinksTotal());
  if (PrintLevel()>2) 
    printf("NIsectLinks=%d\n",fTrackLinkBlock->SiIsectLinkTrk()->NLinksTotal());
  if (fTrackLinkBlock->SiIsectLinkTrk()->fLast+1!=fTrackBlock->NTracks() &&
      fTrackLinkBlock->SiIsectLinkTrk()->fLast!=fTrackBlock->NTracks())
    return 0;
  if (fAreSiHits) {
    if (fTrackLinkBlock->SiHitLinkTrk()->fLast+1!=fTrackBlock->NTracks())
      return 0;
  }
  if (PrintLevel()>2) printf("fLast=%d\n",fTrackLinkBlock->SiHitLinkTrk()->fLast);

  fSiIsectBlock->GetEntry(ientry);
  if (fSiIsectBlock->NSiIsects()<1) return 0;
  if (PrintLevel()>2) printf("NSiIsects=%d\n",fSiIsectBlock->NSiIsects());

  // do whatever you want
  if (PrintLevel()>0) {
    printf(" ---- NSiHits, NTracks: %d, %d\n",nSiHits,fTrackBlock->NTracks());
  }

  TStnLinkBlock* hln=fTrackLinkBlock->SiHitLinkTrk();
  TStnLinkBlock* iln=fTrackLinkBlock->SiIsectLinkTrk();
  if (PrintLevel()>1){
    fSvxDataBlock->Print(0);
    //    fTrackBlock->Print();
    hln->Print(0);
    fSiIsectBlock->Print(0);
    iln->Print(0);
  }

  // Puff up the lookup table for hits in a given digicode
  if (fSvxDataBlock->InitEvent() < 0) {
    printf(" ---- There are no cross references for hits-to-digicode!\n ");
    return 0;
  }

  // Now make hists of hits matched to tracks and wafer intersections
  if (fTreeT || fTreeN || fTreeR) {
    for (int i=0; i<fTrackBlock->NTracks(); i++){
      TStnTrack* t=fTrackBlock->Track(i);
      if (PassTrackCuts(t)) {
	
	if (fAreSiHits) { 
	  if (t->NSvxHits()!=hln->NLinks(i)) {
	    printf("MISMATCH IN TRACKLINKBLOCK!\nNSvxHits=%d, NLinks=%d\n",
		   t->NSvxHits(),hln->NLinks(i));
	    return 0;
	  }
	}

	if (fTreeR) {
	  int b=-1;
	  int prev_b=-2;
	  for (int j=0; j<iln->NLinks(i)*(iln->NLinksTotal()>0); j++){
	    TStnSiIsect* isect = fSiIsectBlock->SiIsect(iln->Index(i,j));
	    if (j==0) prev_b=isect->DigiCode()->Barrel();
	    b=isect->DigiCode()->Barrel();
	    if (b!=prev_b) break;
	  }
	  
	  if (b==prev_b) {
	    fTrack.run=GetHeaderBlock()->RunNumber();
	    fTrack.evt=GetHeaderBlock()->EventNumber();
	    fTrack.b=fSiIsectBlock->SiIsect(iln->Index(i,0))->DigiCode()->Barrel();
	    fTrack.C=2.0*t->C0();
	    fTrack.phi0=t->Phi0();
	    fTrack.d0=t->D0();
	    fTrack.tanl=t->Lam0();
	    fTrack.z0=t->Z0();
	    fTrack.pt=t->Pt();
	    //	    float sd02=(*t->Cov())(3,3);
	    //	    if (sd02<0.0) sd02=9999999.0;
	    //	    fTrack.sd0=std::sqrt(sd02);
	    
	    fHist.fR->Fill();
	  }
	} // make TreeR
	  

	if (fTreeT) {
	  if (hln->NLinksTotal()>0 && t->NSvxHits()>0) {
	    for (int j=0; j<hln->NLinks(i); j++){
	      if (PrintLevel()>1) printf("(i=%d,j=%d)=%d\n",i,j,hln->Index(i,j));
	      TStnSiHit* hit = fSvxDataBlock->SiHit(hln->Index(i,j));
	      int b=hit->DigiCode()->Barrel();
	      int h=hit->DigiCode()->LadderSeg();
	      int p=hit->DigiCode()->PhiWedge();
	      int l=hit->DigiCode()->Layer();
	      int s=hit->DigiCode()->Side();
	      if (PrintLevel()>1){
		printf("(i,j) (b,h,p,l,s) (Q,N): (%d,%d) (%d,%d,%d,%d,%d) (%f,%d)\n",
		       i,j,b,h,p,l,s,hit->Qtot(),hit->Nstrip());
	      }
	      
	      // Fill small ntuple of hits
	      fHit.run=GetHeaderBlock()->RunNumber();
	      fHit.evt=GetHeaderBlock()->EventNumber();
	      fHit.b=b;
	      fHit.h=h;
	      fHit.p=p;
	      fHit.l=l;
	      fHit.s=s;
	      fHit.stp=hit->StripNum();
	      fHit.N=hit->Nstrip();
	      fHit.Q=hit->Qtot()*1.0;
	      fHit.f=std::sin(std::abs(std::atan(1.0/t->Lam0())));
	      fHit.x=hit->Global()->X();
	      fHit.y=hit->Global()->Y();
	      fHit.z=hit->Global()->Z();
	      
	      fHit.phi0=t->Phi0();
	      fHit.z0=t->Z0();
	      fHit.d0=t->D0();
	      fHit.eta=t->Eta();
	      fHit.pt=t->Pt();
	      
	      
	      // Fill the track-hit tree
	      fHist.fT->Fill();
	      
	    } // Loop over SiHit links
	  } // If track has SiHits attached
	} // Make TreeT
	
	if (fTreeN) {
	  for (int j=0; j<iln->NLinks(i)*(iln->NLinksTotal()>0); j++){
	    if (PrintLevel()>1) printf("(i=%d,j=%d)=%d\n",i,j,iln->Index(i,j));
	    TStnSiIsect* isect = fSiIsectBlock->SiIsect(iln->Index(i,j));
	    // Fill small branch of hits nearest to intersection
	    int b=isect->DigiCode()->Barrel();
	    int h=isect->DigiCode()->LadderSeg();
	    int p=isect->DigiCode()->PhiWedge();
	    int l=isect->DigiCode()->Layer();
	    //	int s=isect->DigiCode()->Side();
	    int act=isect->ActiveRegion()->GetBit(0)+
	      isect->ActiveRegion()->GetBit(1)*10+
	      isect->ActiveRegion()->GetBit(2)*100+
	      isect->ActiveRegion()->GetBit(3)*1000;
	    
	    fNear.run=GetHeaderBlock()->RunNumber();
	    fNear.evt=GetHeaderBlock()->EventNumber();
	    fNear.b=b;
	    fNear.h=h;
	    fNear.p=p;
	    fNear.l=l;
	    fNear.act=act;
	    
	    fNear.phi0=t->Phi0();
	    fNear.z0=t->Z0();
	    fNear.d0=t->D0();
	    fNear.eta=t->Eta();
	    fNear.pt=t->Pt();
	    
	    if (fHist.fEff_trk) {
	      if (act==1111) {
		fHist.fEff_trk->Fill(l);
	      }
	    }
	    
	    int smax=1;
	    if (l==0) smax=0;
	    
	    for (int s=0; s<=smax; s++) {
	      isect->DigiCode()->SetSide(s);
	      fNear.s=s;
	      
	      fNear.ix=isect->Global()->X();
	      fNear.iy=isect->Global()->Y();
	      fNear.iz=isect->Global()->Z();
	      if (s==0) {
		fNear.istp=isect->StripNumPhi();
	      }
	      else {
		fNear.istp=isect->StripNumZ();
	      }
	      fNear.ily=isect->LocY();
	      fNear.ilz=isect->LocZ();
	      fNear.f=std::sin(std::abs(std::atan(1.0/t->Lam0())));
	      
	      float d0=7.0E10;
	      float d1=8.0E10;
	      float d2=9.0E10;
	      int m[3];
	      m[0]=-1; m[1]=-1; m[2]=-1;
	      float d;
	      TVector3 x( *(isect->Global()) );
	      int idigi=isect->DigiCode()->fDigiCode;
	      if (fSvxDataBlock->fIndexFirstHit[idigi]>=0 &&
		  fSvxDataBlock->fIndexLastHit[idigi]>=0) {
		for (int k=fSvxDataBlock->fIndexFirstHit[idigi];
		     k<=fSvxDataBlock->fIndexLastHit[idigi]; k++){
		  if (fSvxDataBlock->SiHit(k)->Global()) {
		    d=( *(fSvxDataBlock->SiHit(k)->Global()) - x ).Mag();
		  }
		  else {
		    d=10.E10;
		  }
		  if ( d <= d0 ) {
		    d2=d1;
		    m[2]=m[1];
		    d1=d0;
		    m[1]=m[0];
		    d0=d;
		    m[0]=k;
		  }
		  else if ( d <= d1 ) {
		    d2=d1;
		    m[2]=m[1];
		    d1=d;
		    m[1]=k;
		  }
		  else if ( d <= d2 ) {
		    d2=d;
		    m[2]=k;
		  }
		}
	      }
	      for (int ii=0; ii<3; ii++){
		fNear.stp[ii] = m[ii]>=0 ? fSvxDataBlock->SiHit(m[ii])->StripNum() : -1;
		fNear.N[ii] = m[ii]>=0 ? fSvxDataBlock->SiHit(m[ii])->Nstrip() : -1;
		fNear.Q[ii] = m[ii]>=0 ? fSvxDataBlock->SiHit(m[ii])->Qtot() : -1;
		fNear.x[ii] = m[ii]>=0 ? fSvxDataBlock->SiHit(m[ii])->Global()->X() : 0;
		fNear.y[ii] = m[ii]>=0 ? fSvxDataBlock->SiHit(m[ii])->Global()->Y() : 0;
		fNear.z[ii] = m[ii]>=0 ? fSvxDataBlock->SiHit(m[ii])->Global()->Z() : 0;
	      }

	      // Fill the near-hit tree
	      fHist.fN->Fill();
	    } // Loop over side

	  } // Loop over track intersections
	} // Make TreeN
      } // if desired track
    } // Loop over tracks
  } // make TreeT or TreeN
  if (PrintLevel()>0) {
    printf(" ---- Done with ientry %d\n",ientry);
  }  
  return 0;		       
}

 
//_____________________________________________________________________________
int TSvxAnaModule::EndJob() {
  if (fFile) fFile->Write();
  
  printf("----- end job: ---- %s\n",GetName());
  return 0;
}


bool TSvxAnaModule::PassTrackCuts(TStnTrack* t) {
  bool pass=true;
  pass &= t->Algorithm()==fAlg;
  pass &= t->NSvxHits()>=4;
  pass &= t->NSvxSASHits()>=2;
  pass &= t->NSvxZHits()>=3;

  return pass;
}