////////// ///////// ///////// //////// ///////// ///////// ///////// 72
////////////////////////////////////////////////////////////////////////
// Program name: CDCosmicTracker.cxx
//
// Package: CalDetTracker
//
// Purpose: 
//
// Contact: Chris Smith, Ryan Nichol, Leo Jenner or Jeff Hartnell
////////////////////////////////////////////////////////////////////////

#include "TGraph.h"
#include "TGraphErrors.h"
#include "TROOT.h"
#include "TF1.h"

#include "Plex/PlexSEIdAltL.h"
#include "Plex/PlexStripEndId.h"
#include "MessageService/MsgService.h"

#include "CalDetTracker/CDCosmicTracker.h"
#include "CalDetTracker/CDTrackerOptions.h"

#include <cmath>

using std::cout;
using std::endl;
using std::map;

CVSID("$Id: CDCosmicTracker.cxx,v 1.10 2006/05/27 07:47:41 rhatcher Exp $");

//......................................................................

CDCosmicTracker::CDCosmicTracker()
{
  MSG("CDCosmicTracker",Msg::kVerbose)
    <<"Running default CDCosmicTracker constructor..."<<endl;

  fevenresult=0;
  foddresult=0;

  MSG("CDCosmicTracker",Msg::kVerbose)
    <<"Finished default CDCosmicTracker constructor"<<endl;
}

//......................................................................

CDCosmicTracker::CDCosmicTracker(map<int,CandStripHandle> stripMap,
				 CDTrackerOptions *cdto)
{
  MSG("CDCosmicTracker",Msg::kVerbose)
    <<"Running CDCosmicTracker(Strip) constructor..."<<endl;

  fevenresult = -1; //horizontal strips
  foddresult = -1;  //vertical strips
  fnum_cc_hits = -1;
  fStripMap = stripMap;
  if(!this->SetTrackOptions(cdto)){
    MSG("CDCosmicTracker", Msg::kWarning)
      << "CDTracker Error: No Options Set" << endl;
  }
  MSG("CDCosmicTracker",Msg::kVerbose)
    <<"Finished CDCosmicTracker(Strip) constructor"<<endl;
}

//......................................................................

int CDCosmicTracker::FindStripTrack(int offset)    
{
  MSG("CDCosmicTracker", Msg::kDebug)  
    << "FindTrack( " << offset << ")"<<endl;
  
  if(fStripMap.size()>300) return 2;

  // Get some tracking parameters from the TrackerOptions object  
  int trk_len_cut = fcd_to->GetTrkLenCut();
  int acceptance = fcd_to->GetAcceptance();
  int hits_per_plane = fcd_to->GetHitsPerPlane();
  int max_hot_planes = fcd_to->GetMaxHotPlanes();
  MSG("CDCosmicTracker", Msg::kDebug)  
    << "Track length cut: "<< trk_len_cut << endl;
  
  //temporary map for holding test strips
  map<int,CandStripHandle> pre_tracked;
  //map for holding definite track hits
  map<int,CandStripHandle> tracked; 
  
  StripEnd::EStripEnd dum_strip_end_crate = StripEnd::kWhole;

  //Set up some large arrays to make TGraphs
  float planeHits[500] = {};
  float stripHits[500] = {};
  float planeErrs[500] = {};
  float stripErrs[500] = {};
  
  float planeHits2[500] = {};
  float stripHits2[500] = {};
  float planeErrs2[500] = {};
  float stripErrs2[500] = {};

  //Some tracking variables
  int numHitsInSection2=0;
  int hitNum=0;
  int lastPlaneHit=-1;
  int numHitsInSection=0;

  //Fill arrays:
  for(int i=0;i<60;i++){
    for(int j=0;j<24;j++){

      if(i%2==offset) {
	PlexStripEndId dummyplex(Detector::kCalDet,
				 i,j,dum_strip_end_crate);
	int theKey = dummyplex.Build18BitPlnStripKey();	
	if(fStripMap.find(theKey)!=fStripMap.end()){
	  
	  MSG("CDCosmicTracker", Msg::kDebug)  
	    << "Really Before: "<< i << " " << j << endl;
	  
	  if(lastPlaneHit>-1 && i-lastPlaneHit>6) {
	    MSG("CDCosmicTracker", Msg::kDebug)  
	      << "New Section" << endl;
	    
	    //This is a new section.
	    if(numHitsInSection>2*(trk_len_cut-1)) {
	      MSG("CDCosmicTracker", Msg::kDebug)  
		<< "With: " << numHitsInSection << " or " << hitNum 
		<< " hits" << endl;
	      
	      //Then we will consider it a valid section for now.
	      for(int index=0;index<hitNum;index++) {
		planeHits2[index]=planeHits[index];
		planeErrs2[index]=planeErrs[index];
		stripHits2[index]=stripHits[index];
		stripErrs2[index]=stripErrs[index];
		
		MSG("CDCosmicTracker", Msg::kDebug)  
		  << "In between: " << planeHits2[index] << " "
		  << stripHits2[index] << endl;
		
		planeHits[index]=0;
		planeErrs[index]=0;
		stripHits[index]=0;
		stripErrs[index]=0;
	      }
	      numHitsInSection2=numHitsInSection;
	      numHitsInSection=0;
	      hitNum=0;
	    }
	    else {
	      for(int index=0;index<hitNum;index++) {
		planeHits[index]=0;
		planeErrs[index]=0;
		stripHits[index]=0;
		stripErrs[index]=0;
	      } 
	      numHitsInSection=0;
	      hitNum=0;
	    }
	  }
	  


	  if(fStripMap[theKey].GetCharge(CalDigitType::kPE,
					 StripEnd::kPositive)>0){
	    planeHits[hitNum]=i;  // z position in units of "plane"
	    planeErrs[hitNum]=0.; // set error on z position to be zero
	    stripHits[hitNum]=j;  // transverse position in units of "strip"
	    //** error on transverse position set to be 1./sqrt(adc)
	    stripErrs[hitNum]= 1.00/
	      sqrt(fStripMap[theKey].GetCharge(CalDigitType::kPE,
					       StripEnd::kPositive));
	    
	    hitNum++;
	    numHitsInSection++;
	    lastPlaneHit=i;
	  }
	  
	  if(fStripMap[theKey].GetCharge(CalDigitType::kPE,
					 StripEnd::kNegative)>0){
	    planeHits[hitNum]=i;  // z position in units of "plane"
	    planeErrs[hitNum]=0.; // set error on z position to be zero
	    stripHits[hitNum]=j;  // transverse position in units of "strip"
	    //** error on transverse position set to be 1./sqrt(adc)
	    stripErrs[hitNum]= 1.00/
	      sqrt(fStripMap[theKey].GetCharge(CalDigitType::kPE,
					       StripEnd::kNegative));
	    
	    hitNum++;
	    numHitsInSection++;	
	    lastPlaneHit=i;
	  }
	}
      }
    }
  }
  
  int whichPlanes[60]={0};
  int whichPlanes2[60]={0};
  for(int index=0;index<numHitsInSection;index++) {
    whichPlanes[int(planeHits[index])]++;
  }
  
  for(int index=0;index<numHitsInSection2;index++) {
    whichPlanes2[int(planeHits2[index])]++;
  }

  int numPlanes=0;
  int numPlanes2=0;
  for(int index=0;index<60;index++) {
    if(whichPlanes[index]) {
      numPlanes++;
      MSG("CDCosmicTracker", Msg::kDebug) 
	<<"Got Plane: "<<index<<" with "<<whichPlanes[index]<<endl;
    }
    if(whichPlanes2[index]) {
      numPlanes2++;
      MSG("CDCosmicTracker", Msg::kDebug) 
	<<"Got Plane2: "<<index<<" with "<<whichPlanes2[index] << endl;
    }
  }
  
  if(numPlanes2>numPlanes) { 
    MSG("CDCosmicTracker", Msg::kDebug)  
      <<"There are: " << numHitsInSection << " hits & " 
      << numHitsInSection2 << " other hits" <<  endl;
    
    //Simple definition but will do for now.
    for(int index=0;index<numHitsInSection2;index++) {
      MSG("CDCosmicTracker", Msg::kDebug)  
	<< "Middle: " << planeHits2[index] << " " 
	<< stripHits2[index] << endl;
      
      planeHits[index]=planeHits2[index];
      planeErrs[index]=planeErrs2[index];
      stripHits[index]=stripHits2[index];
      stripErrs[index]=stripErrs2[index];
    }
    numHitsInSection=numHitsInSection2;
    hitNum=numHitsInSection2;  
  }
  
  if(hitNum<=3) return 0;
  
  TGraphErrors gr(hitNum,planeHits,stripHits,planeErrs,stripErrs);
  TF1 fitty("fitty","pol1",-1,60);
  gROOT->SetBatch(1);
  gr.Fit("fitty","QR");
  gROOT->SetBatch(0);
  
  int tempLastHits=hitNum;
  int newTempNumHits=0;
  float newPlaneHits[500]={};
  float newStripHits[500]={};
  float newPlaneErrs[500]={};
  float newStripErrs[500]={};
  float tempPlaneHits[500]={};
  float tempStripHits[500]={};
  float tempPlaneErrs[500]={};
  float tempStripErrs[500]={};
  
  for(int i=0;i<tempLastHits;i++) {
    MSG("CDCosmicTracker", Msg::kDebug)  << "Before: " << planeHits[i] << " " 
					 << stripHits[i] << endl;
    tempPlaneHits[i]=planeHits[i];
    tempPlaneErrs[i]=planeErrs[i];
    tempStripHits[i]=stripHits[i];
    tempStripErrs[i]=stripErrs[i];
  }
  
  while(newTempNumHits!=tempLastHits) {
    newTempNumHits=0;     
    
    for(int i=0;i<tempLastHits;i++) {
      float diffFromLine=fabs(fitty.Eval(tempPlaneHits[i])
			      -tempStripHits[i]);
      if(diffFromLine<3) {
	newPlaneHits[newTempNumHits]=tempPlaneHits[i];
	newPlaneErrs[newTempNumHits]=tempPlaneErrs[i];
	newStripHits[newTempNumHits]=tempStripHits[i];
	newStripErrs[newTempNumHits]=tempStripErrs[i];
	newTempNumHits++;
      }
    }
    
    if(newTempNumHits<=3) break;
    
    tempLastHits=newTempNumHits;
    
    for(int i=0;i<tempLastHits;i++) {
      tempPlaneHits[i]=newPlaneHits[i];
      tempPlaneErrs[i]=newPlaneErrs[i];
      tempStripHits[i]=newStripHits[i];
      tempStripErrs[i]=newStripErrs[i];
    }
    
    TGraphErrors gr2(newTempNumHits,newPlaneHits,
		     newStripHits, newPlaneErrs,
		     newStripErrs);
    gROOT->SetBatch(1);
    gr2.Fit("fitty","QR");
    gROOT->SetBatch(0);
    
  }
  
  if(newTempNumHits<=3) return 0;
    
  int piShowerFlag=0;  // if plane has more than 2 hits in tracking 
                       // window increment pi_shower_flag
  int isAMuon=1;
  int firstHit=-1;
  int lastHit=-1;
  
  int numTrackedHits=0;
  
  for(int i=0;i<60;i++){
    
    int planeNumHits=0;
    int gotThisPlane=0;
    
    for(int index=0;index<hitNum;index++) {

      if(i==planeHits[index]) {
	gotThisPlane=1;
	break;
      }
      
    }
    
    if(!gotThisPlane) continue;
    
    for(int j=0;j<24;j++){
      
      if(i%2==offset) {
	
	PlexStripEndId dummyplex3(Detector::kCalDet,i,
				    j,dum_strip_end_crate);
	int theKey3 = dummyplex3.Build18BitPlnStripKey();	
	
	if(fStripMap.find(theKey3)!=fStripMap.end()){
	  
	  if(firstHit==-1) firstHit=i;
	  
	  if((j>fitty.Eval(i)-acceptance) && (j<fitty.Eval(i)+acceptance))
	    planeNumHits++;
	  
	  if((j>fitty.Eval(i)-1.5) && (j<fitty.Eval(i)+1.5)) {
	    if(offset==0) fEvenPlStrips[theKey3] = fStripMap[theKey3];
	    else fOddPlStrips[theKey3] = fStripMap[theKey3];
	    lastHit=i;
	    if(fStripMap[theKey3].
	       GetCharge(CalDigitType::kPE,StripEnd::kNegative)>0) 
	      numTrackedHits++;
	    if(fStripMap[theKey3].
	       GetCharge(CalDigitType::kPE,StripEnd::kPositive)>0) 
	      numTrackedHits++;
	  }  
	}
      }
    }
    
    if(planeNumHits >= (hits_per_plane*2)) isAMuon=2;
    else if(planeNumHits>4) piShowerFlag++;
    else piShowerFlag = 0;
    
    if(piShowerFlag>max_hot_planes) isAMuon=2;   
  }
  
  if(lastHit-firstHit<(2*(trk_len_cut-1))) isAMuon=0;
  if(hitNum>1.7*numTrackedHits) isAMuon=4;
  
  MSG("CDCosmicTracker", Msg::kDebug)  
    << "Last and First for offset: " 
    << offset << " " << lastHit 
    <<" "<< firstHit <<endl;
  return isAMuon;

}

//......................................................................