/*!
 * \class StGeantEleFinderMaker
 * \author Naresh Subba KSU,  Jan Balewski IUCF , October 2006
 * Goal: identify electrons from the Geant record
 * $Id: $
 */


#include <TH2.h>
#include <StMessMgr.h>

#include "StGeantEleFinderMaker.h"

#include "StChain.h"
#include "St_DataSetIter.h"
#include "StMcEventMaker/StMcEventMaker.h"
#include "StMcEvent/StMcEvent.hh"
#include "StMcEvent/StMcVertex.hh"
#include "StMcEvent/StMcTrack.hh"
#include "StMcEvent/StMcTpcHit.hh"
#include "StMcEvent/StMcFstHit.hh"
#include "StMcEvent/StMcIgtHit.hh"
//#include "StMcEvent/StMcHit.hh"

ClassImp(StGeantEleFinderMaker)

  struct JanHitZ_LT {
    bool operator() ( StMcHit* lhs,StMcHit* rhs) {return lhs->position().z() < rhs->position().z();}
  };   


StGeantEleFinderMaker::StGeantEleFinderMaker(const char *name):StMaker(name){
  memset(hA,0,sizeof(hA));
   
  HList=0;
  par_pTmin=5.;
 
  gMessMgr->Message("","I")
    << "GeantEleFinder::cuts\n"
    << Form("pT>%.2f GeC/c\n",par_pTmin) 
    << endm;
  
}


StGeantEleFinderMaker::~StGeantEleFinderMaker(){
  //nop
}


//_____________________________________________________________________________
//
Int_t 
StGeantEleFinderMaker::Init(){
  
  //.................. histos ....................
  // histo naming convention:
  int nEtaB=34;
  float eta1=0.6, eta2=2.7;

  hA[0]=new TH1F("gEveStat","Geant Event counter",30,0.5,30.5);
  hA[1]=new TH1F("gTrkStat","Geant track counter",20,0.5,20.5); 
  
  hA[2]=new TH2F("gZ_Eta_fst", "Geant hits in FST; gen #eta; Z hit (cm)",nEtaB,eta1,eta2,75,25,40);
  hA[3]=new TH2F("gNfst","# of FST hits ;  gen #eta; N hit",nEtaB,eta1,eta2,6,-0.5,5.5);

  hA[4]=new TH2F("gZ_Eta_igt", "Geant hits in IGT; gen #eta; Z hit (cm)",nEtaB,eta1,eta2,55,50,160);
  hA[5]=new TH2F("gNigt","# of IGT hits ;  gen #eta; N hit",nEtaB,eta1,eta2,6,-0.5,5.5);

  hA[6]=new TH2F("gR_Eta_tpc", "Geant hits in TPC; gen #eta; Rxy hit (cm)",nEtaB,eta1,eta2,150,50,200);
  hA[7]=new TH2F("gNtpc","# of TPC hits ;  gen #eta; N hit",nEtaB,eta1,eta2,50,-0.5,49.5);


  hA[8]=0;
  hA[9]=new TH2F("gNall","# of VER+FST+IGT+TPC hits ;  gen #eta; N hit",nEtaB,eta1,eta2,60,-0.5,59.5);


  hA[10]=new TH1F("gPtG", "Geant  prim tracks; pT(GeV/c) ",50,0.,10.);
  hA[11]=new TH1F("gEtaG", "Geant  prim  tracks; #eta ",40,0.,2.8); 
  hA[12]=new TH1F("gPhiG", "Geant  prim  tracks; #phi (rad)",60,-4.,4.);
  hA[13]=new TH1F("gEneG", "Geant  prim tracks; energy (GeV) ",50,0.,20.);
  

  hA[20]=new TH2F("gL1_Eta_all", "g track L1=len(lastH-firstH) , SFT+IGT+TPC; gen #eta; track len (cm)",nEtaB,eta1,eta2,50,0,220);

  hA[21]=new TH2F("gL2_Eta_all", "g track L2=len(vert - midH) , SFT+IGT+TPC; gen #eta; track len (cm)",nEtaB,eta1,eta2,50,0,220);
  hA[22]=new TH2F("gVf_Eta_all", "g vertex factor  Vf=(L2-L1/2)/L1; gen #eta; vertex factor",nEtaB,eta1,eta2,50,0,2.);
  
  //... add all histos to the list
  assert(HList);
  int i;
  for(i=0;i<mxHa;i++) if(hA[i]) HList->Add(hA[i]);
  
  return StMaker::Init();
}

//________________________________________________
//________________________________________________
void
StGeantEleFinderMaker::Clear(const Option_t*){
  hitL.clear();
}

//________________________________________________
//________________________________________________
Int_t 
StGeantEleFinderMaker::Make(){
  hA[0]->Fill(1);
  float gEta;
  findMcTrk(gEta);
  if(hitL.size()>0) hA[0]->Fill(2);
  
  printf("GeantEleFinder hitL size=%d\n",hitL.size());
  uint i;

  Hit3D *h1=0, *h2=0, *hv=0;
  for(i=0;i<hitL.size();i++) {
    Hit3D *h=&hitL[i];
    printf("i=%d type=%d  x,y,z=%.2f %.2f %.2f\n",i, h->type,   h->r.x(),   h->r.y(),   h->r.z());
    if(h->type==Hit3D::VER) hv=h;
    
    if(h->type==Hit3D::VER) continue;
    //.... only hits in the trackers, find first & last
    if(h1==0) h1=h;
    if(h2==0) h2=h;
    if(h1->r.z() >  h->r.z() ) h1=h; // left end
    if(h2->r.z() <  h->r.z() ) h2=h; // right end
  }
  if(h1==0) return kStOK;  // no hits in  trackers
  if(h1==h2) return kStOK;  // one hit in  trackers

  printf("Lhit type=%d  x,y,z=%.2f %.2f %.2f\n", h1->type,   h1->r.x(),   h1->r.y(),   h1->r.z());
  printf("Rhit type=%d  x,y,z=%.2f %.2f %.2f\n", h2->type,   h2->r.x(),   h2->r.y(),   h2->r.z());

  TVector3 diff1=h2->r - h1->r;// approximated track length (assumes traight line)
  float Len1=diff1.Mag(); 
  printf("Len1=%.2f\n",Len1);
  hA[20]->Fill(gEta,Len1);

  // vertex properties
  TVector3 diff2=0.5*(h1->r+ h2->r) - hv->r; // from vertex to middel of tracker track
  float Len2=diff2.Mag(); 
  printf("Len2=%.2f\n",Len2);
  hA[21]->Fill(gEta,Len2);

  float vertFac=-999;
  if(Len1>0)  vertFac=(Len2-Len1/2.)/Len1;
  hA[22]->Fill(gEta,vertFac);
  printf("vertFac=%.3f\n",vertFac);


  return kStOK;  
}


//________________________________________________
//________________________________________________
void
StGeantEleFinderMaker::findMcTrk(float &gEta){  
  // **************************************
  // get pointer to mcEventMaker, Event *
  // **************************************

  float par_fst_zStep=0.5; //cm
  float par_igt_zStep=0.5; //cm
  float par_tpc_zMax=198; //cm
  gEta=-999;
  StMcEvent* mMcEvent = 0;  //!
  mMcEvent = (StMcEvent*) StMaker::GetChain()->GetDataSet("StMcEvent"); 
  assert(mMcEvent); 
  // primary vertices
  const StSPtrVecMcVertex &VL= mMcEvent->vertices();
  
  printf("%s::nVert=%d, scan only the first\n", GetName(), VL.size());
  uint i;
  for(i=0;i<VL.size();i++) { 
    if(i>0) break; // ignore 2ndary vertices
    StMcVertex* V=VL[i];
    // Rxy vx z vertex to remove fake electrons
    float x=V->position().x();
    float y=V->position().y();
    float z=V->position().z();
    float Rxy=sqrt(x*x+y*y);
    int nD=V-> numberOfDaughters();    
    printf("vx=%f,vy=%f,vz=%f Rxy=%f  #daughters=%d\n",x,y,z,Rxy,nD);

    Hit3D hit;
    hit.r.SetXYZ(V->position().x(), V->position().y(), V->position().z());
    hit.type=Hit3D::VER;
    hitL.push_back(hit);  
    
    for(int it=0; it<nD;it++ ) { // primary tracks
      assert(it==0); // should be 1 prim electron per event
      StMcTrack* tr=V->daughter(it);
      float phi=tr->momentum().phi(); 
      gEta=tr->pseudoRapidity();              
      float pt=tr->pt();     
      int gid=tr->geantId();
      hA[1]->Fill(1); 
      if(pt<par_pTmin) continue;
      hA[1]->Fill(2); 
      //if(gid!=3 && gid!=2) continue;
      assert(gid==3);
      // if(eta>1.4) continue;//tmp

      int nFstH=tr->fstHits().size();  // geant FST hits
      int nIgtH=tr->igtHits().size();  // geant IGT hits
      int nTpcH=tr->tpcHits().size();  // geant TPC hits
      int naFstH=0,naIgtH=0,naTpcH=0; // # of accepted hits 


      printf("g2t: i=%d pt=%f phi=%f eta=%f nFstH=%d nIgtH=%d  nTpcH=%d \n",it,pt,phi,gEta,nFstH,nIgtH,nTpcH);

      float zLast=-999;

      //........... FST........................
      std::sort(tr->fstHits().begin(), tr->fstHits().end(), JanHitZ_LT()); // sort hits vs. Z

      for(int ih=0; ih<nFstH;ih++) {
	StMcFstHit *h=tr->fstHits()[ih];
	float zNew=h->position().z();
	assert(zNew>=zLast); // tmp
	//	printf("i=%d %f %f\n",ih,zLast,zNew);
	if(zLast+par_fst_zStep > zNew) continue;
	zLast=zNew;
	hA[2]->Fill(gEta,zNew);
	naFstH++;
	
	Hit3D hit;
	hit.r.SetXYZ(h->position().x(), h->position().y(), h->position().z());
	hit.type=Hit3D::FST;
	hitL.push_back(hit);  
	//	printf("McFstHit lay=%d ladd=%d pln=%d mod=%d x,y,z=%.2f %.2f %.2f\n",(int) h->layer(),(int) h->ladder(),(int) h->plane(),(int) h->module(), h->position().x(), h->position().y(), h->position().z());

      }
      
      hA[3]->Fill(gEta,naFstH);
      
      //........... IGT........................
      std::sort(tr->igtHits().begin(), tr->igtHits().end(), JanHitZ_LT()); // sort hits vs. Z
 
      zLast=-999; // reset it just in case
      for(int ih=0; ih<nIgtH;ih++) {
	StMcIgtHit *h=tr->igtHits()[ih];
	float zNew=h->position().z();
	//	printf("i=%d %f %f\n",ih,zLast,zNew);
	if(zLast+par_igt_zStep > zNew) continue;
	zLast=zNew;
	hA[4]->Fill(gEta,h->position().z());
	naIgtH++;
	assert(zNew>=zLast); // tmp
	Hit3D hit;
	hit.r.SetXYZ(h->position().x(), h->position().y(),h->position().z());
	hit.type=Hit3D::IGT;
	hitL.push_back(hit);  

	//	if(naIgtH==1) printf("\n zLast=%f\n",zLast);
	//	printf("McIgtHit lay=%d ladd=%d waf=%d side=%d x,y,z=%.2f %.2f %.2f\n",(int) h->layer(),(int) h->ladder(),(int) h->wafer(),(int) h->side(), h->position().x(), h->position().y(), h->position().z());	
      }
      hA[5]->Fill(gEta,naIgtH);

  
      //........... TPC........................
      float secLast=-1, padLast=-1;
      for(int ih=0; ih<nTpcH;ih++) {
	StMcTpcHit *h=tr->tpcHits()[ih];
	float zNew=h->position().z();
	int secNew=h->sector();
	int padNew=h->padrow();
	if(padNew==13) continue; // dead padrow 13
	if(padNew==1) continue; //  padrow 13

	if( fabs(zNew)>par_tpc_zMax) continue; // clip large z	
	if(secNew==secLast && padNew==padLast) continue; // drop double hits
	secLast=secNew;
	padLast=padNew;

	naTpcH++;
	Hit3D hit;
	hit.r.SetXYZ(h->position().x(), h->position().y(),h->position().z());
	hit.type=Hit3D::TPC;
	hitL.push_back(hit);  

	float x= h->position().x();
	float y= h->position().y();
	float phi=(3-secNew)*3.1416/6.;
	float rN=x*cos(phi)+y*sin(phi); // normal direction
	hA[6]->Fill(gEta,rN); // WARN Rxy is defined somewhere?! 
	
	//if(naTpcH==1) printf("\n");
	//printf("McTpcHit sec=%d padrow=%d  x,y,z=%.2f %.2f %.2f  rN=%.2f\n",secNew, padNew , h->position().x(), h->position().y(), h->position().z(),rN);
      }
     hA[7]->Fill(gEta,naTpcH);

      
      printf("g2t: accepted i=%d pt=%f phi=%f eta=%f nFstH=%d nIgtH=%d  nTpcH=%d \n",it,pt,phi,gEta,naFstH,naIgtH,naTpcH);

      hA[9]->Fill(gEta,hitL.size());
      
      // fill all tracks,no cut yet


      // all tracks are accepted from here ....
      
      // add track to the list
      

      
      //....  QA histos
      hA[10]->Fill(pt);
      hA[11]->Fill(gEta);
      hA[12]->Fill(phi);  
      hA[13]->Fill(tr->energy());  
      
    } //end of loop over tracks in a single vertex
  } //end of loop over vertices
}



// $Log: $