// Run Fast MC interactively with Root
// Sample program using TrjPart              Sep/23/99    M. Iwasaki
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x Trk_trj.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x Trk_trj.C(300)
//

#include <stdlib.h>

int Trk_trj(int   nEvent = 10)
{ 
  gROOT->Reset();

  // In order to refer to shareable libraries in this simple form,
  // include their paths in a .rootrc file.
  gSystem->Load("libLCDEvent.so");
  gSystem->Load("libLCDRootAppsUtil.so");
  gSystem->Load("libLCDFastMC.so");
  gSystem->Load("libLCDPhUtil.so");

  // Input HEPEVT file
  char* stdFile = "test.dat";

  // Detector parameter file
  char* geomFile = 
    "/afs/slac/g/nld/RootApps/V2.1/ParFiles/Small.par";

  // Tracker smearing file
  char* smearFile = 
    "/afs/slac/g/nld/RootApps/V2.1/ParFiles/lookup_small.nobmcon";

  LCDreadStdFile *source = new LCDreadStdFile(stdFile);
  LCDFastMC* fmc = new LCDFastMC(source, geomFile, smearFile, 0);


  // Define a histogram or two
  TH1F *dL = new TH1F("dL","Minimum clus-trk distance",100,0.,20.0);
  dL->SetFillColor(4);

  // Set up TrajPart!
  LCDTrjPart trjpart;
  LCDTrjPart* ptrjb=&trjpart;
  ptrjb->GetSmall();

  const Double_t mpi     = 0.1396;

  // Event loop
  for (Int_t iEvent = 0; iEvent < nEvent; iEvent++) {

    LCDEvent* event = new LCDEvent();
    event->Create();

    Int_t ierr = fmc->FetchEvent(event);
    if (!ierr) break;

    fmc->Cleanup();
    fmc->Doit(event);
    
    // McPart
    TObjArray* McPart = event->MCparticles();
    Int_t nMCpart = McPart->GetEntries();

    

    //charged tracks
    TObjArray* Tracks = event->Tracks();
    Int_t nTRKS = Tracks->GetEntries();

    Double_t Etrk, Ptrk, COStrk;
    Double_t cal_pos[3] ;
    Double_t momentum[4];
    Double_t position[3];
    
    //Neutral clusters
    TObjArray* Clusters = event->ClusterLst();
    Int_t nCluster = Clusters->GetEntries();

    for (Int_t iClus = 0; iClus < nCluster; iClus++) {
      LCDCluster *Clus = (LCDCluster *) Clusters->At(iClus);
      if (Clus) {

	Double_t     e     =  Clus->GetEnergy();
	Double_t     Phi   =  Clus->GetEnergyPhi();
	Double_t     Theta =  Clus->GetEnergyTheta();
	Double_t     R_cal =  Clus->GetEnergyR();

	// The definition r = distance from IP 
	R_cal = R_cal*sin(Theta);
	
        cal_pos[0] = R_cal * cos(Phi);
	cal_pos[1] = R_cal * sin(Phi);
	cal_pos[2] = R_cal * cos(Theta)/sin(Theta);

	// Cluster cut
        if (e < 0.1 || fabs(cos(Theta))> 0.90 ) continue;

	Int_t barend = -1 ; Double_t rt;
	
        // Barrel or Endcap :: Small
	if((R_cal> 75.&&R_cal<140.&& fabs(cal_pos[2])<185.) || //EM  Barrel 
	   (R_cal>140. ))                                   { //HAD Barrel 
	    barend = 0;
	    rt = R_cal;
        }
	else if( (R_cal< 75.0 && fabs(cal_pos[2])<185.0) || //EM  endcap
	    (R_cal<140.0 && fabs(cal_pos[2])>185.0) ){ //HAD endcap
	    barend = 1;
	    rt = cal_pos[2];
	} 

	if ( barend == -1 ) {
	  printf("barend = -1  r, z = %f %f \n", R_cal, cal_pos[2]);
	  continue;
	}

	//cluster -track matching 
	Double_t d_L_min = 1000000.; 
	for (Int_t iTRK = 0; iTRK < nTRKS; iTRK++) {
          LCDTrack* TRK = (LCDTrack *) Tracks->At(iTRK);
          if (TRK) {
	    Double_t* TRK_mom = TRK->GetMomentum();	
	    Ptrk   = 0.;
	    for(Int_t i=0; i<3; i++ ) { 
	      Ptrk += (*(TRK_mom+i))*(*(TRK_mom+i)); }
	    Etrk   = Ptrk + mpi* mpi;
	    Etrk   = sqrt(Etrk); Ptrk   = sqrt(Ptrk);
	    COStrk = (*(TRK_mom+2))/Ptrk; 
	    // Track cut
	    if ( Etrk < 0.2 ||  fabs(COStrk) > 0.90) continue;

	    Double_t* TRK_pos = TRK->GetPosition();	
	    Int_t    trk_charge = TRK->GetCharge();

	    // Get the extrapolated point
	    for (Int_t i=0; i<3;i++){
	      momentum[i] =  (*(TRK_mom+i));
	      position[i] =  (*(TRK_pos+i));
	    }
	    momentum[3] = Etrk;

	    Double_t Proj = cal_pos[0]*momentum[0] +
	      cal_pos[1]*momentum[1] +cal_pos[2]*momentum[2];	  
	    if (Proj < 0.) continue;
	    
//  Trajection
	    ptrjb->SetTrack(McPart,trk_charge,momentum,position,barend,rt);    

	    Double_t*  xtrj=ptrjb->GetPos();
	    
	    Double_t d_L=0.;
	    Double_t d_theta;
	    if (barend == 0){  // Barrel
		d_theta = xtrj[0]*cal_pos[0] + xtrj[1]*cal_pos[1];
		d_theta = d_theta / sqrt(xtrj[0]*xtrj[0] + xtrj[1]*xtrj[1]);
		d_theta = d_theta / sqrt(cal_pos[0]*cal_pos[0] + 
					 cal_pos[1]*cal_pos[1]);

		//		printf("d_theta = %f %f %f \n",
		//       d_theta, xtrj[0], xtrj[1]);

		if (d_theta < 1. && d_theta > -1.) 
		  { d_theta = acos(d_theta);}
		else if (d_theta >=  1.)
		  {printf("dcos >=  1,dcos = %f \n",d_theta); d_theta = 0.;} 
		else 
		  {printf("dcos <= -1,dcos = %f \n",d_theta); d_theta = 3.142;}
		
		d_L = sqrt((rt*d_theta)*(rt*d_theta) + 
			   (xtrj[2]-cal_pos[2])*(xtrj[2]-cal_pos[2]));
	    }
	    else if (barend == 1){ // EndCap
		d_L = 0.;
		for (Int_t ii=0;ii<2;ii++){
		    d_L+=(xtrj[ii]-cal_pos[ii])*(xtrj[ii]-cal_pos[ii]);
		}
		d_L = sqrt(d_L);
		if ( cal_pos[2]*(*(xtrj+2)) < 0) d_L = 99999.; 	     
	    }
	    
	    if (d_L < d_L_min) d_L_min = d_L ; 
     
	  } // if (TRK) 
	} // track loop
        
	// Minimum trk-cluster distance
	dL->Fill(d_L_min);
	
      } //if (Clus)
    }//cluster loop

    event->Clean();    delete event;
    event = 0;

  } // event loop
  

  TCanvas *c1 = new TCanvas("c1", "c1", 400,400);
  c1->Divide(1, 1);
  c1->cd(1);  dL->Draw("B");
  return 1;
}