// ----------------------------------------------------------------------------
// $Id: stdhep_GoTrkTrj.C,v 1.6 2001/09/20 01:47:53 toshi Exp $
// ----------------------------------------------------------------------------
//
// Run Fast MC interactively with Root
// Sample program using TrjPart              04/26/2001    Toshi
//
// :: How to use ::
// 
// 1) type root
// 2)
//     root [0]  .x stdhep_GoTrkTrj.C
// 
// .. that's it.
//
// For example, if you like to generate 300 events, type like
//     root [0]  .x stdhep_GoTrkTrj.C(300)
//


Int_t stdhep_GoTrkTrj(Int_t nEvent = 100) {
  gROOT->Reset();

  TString parfile_dir;
  parfile_dir += gSystem->Getenv("LCDROOT")   ;parfile_dir += "/";
  parfile_dir += "ParFiles/";

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

  LCDEvent event;

  // Input generated (stdHEP format) file
  Char_t* stdFile = "test.dat";
  LCDreadStdFile source(stdFile,&event);

  // Open detector parameter file
  Char_t* geomFile = 
    "SDMar01.par"; // SD
  //"LDMar01.par"; // LD
  TString s_geomFile(parfile_dir);
  s_geomFile += geomFile;
  FILE* m_geomfile =  fopen(s_geomFile.Data(),"r");
  LCDGetParameters gp(m_geomfile);
  fclose(m_geomfile);

  // If you want to change some detector parameter(s) for FastMC 
  // do here
  // For example..
  //gp.SetMagneticField(2.);  // Change Magnetif field
  //gp.SetCoilinHAD();        // Set HAD inside the Coil (This case
                              // you need to change some detector size also..)

  // Tracking smearing parameter file
  Char_t* smearFile = 
    "lookup_silicon.nobmcon";
  //"lookup_large.nobmcon";
  TString s_smearFile(parfile_dir);
  s_smearFile += smearFile;

  // Initialize FastMC
  LCDFastMC fmc(&gp,s_smearFile.Data(),&event); 
  // If you don't want the Cluster merging
  //fmc.SetNoClusterMerge();

  // 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 ptrj(&gp);

  const Double_t mpi     = 0.1396;

  // Event loop
  Int_t iEvent;
  for (iEvent = 0; iEvent < nEvent; iEvent++) {
    if (!source.GetEvent()) break; // Read an event from the stdhep file.
    fmc.Doit();                    // Through FastMC
    
    //charged tracks
    TObjArray* Tracks = event->Tracks();
    Int_t nTRKS = Tracks->GetEntries();

    Double_t Etrk, Ptrk, COStrk;
    TVector3 cal_pos(0.,0.,0.) ;
    TVector3 momentum(0.,0.,0.);
    TVector3 position(0.,0.,0.);
    TVector3 xtrj(0.,0.,0.);

    //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();
	Int_t        barend=  Clus->GetBarEnd();

	// The definition r = distance from IP 
	R_cal = R_cal*sin(Theta);
	
        cal_pos.SetX( R_cal * cos(Phi) );
	cal_pos.SetY( R_cal * sin(Phi) );
	cal_pos.SetZ( R_cal * cos(Theta)/sin(Theta) );

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

	Double_t rt;	
	if(barend == 0) { rt = R_cal; }      // Barrel
        else            { rt = cal_pos.Z(); }// EndCap

	//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) {
	    momentum = TRK->GetMomentumVector(0.);	
	    Ptrk   = momentum.Mag();
	    Etrk   = TMath::Sqrt(TMath::Power(Ptrk,2) + 
				 TMath::Power(mpi,2)   );
	    COStrk = momentum.CosTheta();
	    // Track cut
	    if ( Etrk < 0.2 ||  TMath::Abs(COStrk) > 0.90) continue;

	    position = TRK->GetPositionVector(0.);	
	    Int_t    trk_charge = TRK->GetCharge();

	    Double_t Proj = cal_pos * momentum;
	    if (Proj < 0.) continue;
	    
//  Trajection
	    ptrj.Swim(trk_charge,momentum,position,barend,rt);    

	    xtrj=ptrj.GetPos();
	    
	    Double_t d_L=0.;
	    Double_t d_theta, d_theta2;
	    if (barend == 0){  // Barrel
	      d_theta = xtrj.X()*cal_pos.X() + xtrj.Y()*cal_pos.Y();
	      d_theta = d_theta / xtrj.Pt();
	      d_theta = d_theta / cal_pos.Pt();
	      
	      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 = TMath::Pi();
	      }
	      
	      d_L = sqrt((rt*d_theta)*(rt*d_theta) + 
			 (xtrj.Z()-cal_pos.Z())*(xtrj.Z()-cal_pos.Z()));
	    }
	    else if (barend == 1){ // EndCap
	      d_L = (xtrj - cal_pos).Mag();
	      if ( cal_pos.Z() * xtrj.Z() < 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 loop
  

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