////////////////////////////////////////////////////////////////////////////////////
//Macro to make a root file from the info of tracks contained in an IcePack file
//LB IIHE,VUB Brussels 23/10/2013
///////////////////////////////////////////////////////////////////////////////////

#include "/user/lbrayeur/codes/Analysis/Cuts/src/utils_include.h"
#include "/user/lbrayeur/codes/Analysis/Cuts/src/Cuts.h"

using namespace std;

int main(int argc, const char* argv[]){

  //Load data
  TString filename = argv[1] ;

  TChain* data=new TChain("T");
 
  data->Add(filename);
  
  //Define the type of data
  TString argtype=argv[2];
  Int_t type=0;
  if(argtype=="cors")type=1;
  if(argtype=="data")type=2;

  //Set the output file
  TString OUTfilename = argv[3] ;

  //The different reco
  Int_t NrecoHMB=3;
  Int_t NrecoB=3;
  Int_t NrecoRL=3;
  TString* recoHMB=new TString[NrecoHMB];
  TString* recoB=new TString[NrecoB];
  TString* recoRL=new TString[NrecoRL];
  recoHMB[0]="IceDwalk";
  recoHMB[1]="IcePandel4Dwalk";
  recoB[0]="IceLinefit4Dwalk";
  recoB[1]="IceLinefitI";
  recoHMB[2]="IcePandel4LinefitI";
  recoRL[0]="SPEFit2";
  recoRL[1]="MPEFit";
  recoRL[2]="MuonLlhFit";
  recoB[2]="LineFit";
  
  
 
  // Define a pointer for an event
  IceEvent* evt=new IceEvent();

  // Branch in the tree for the event input
  data->SetBranchAddress("IceEvent",&evt);


  cout << endl;
  cout << " *READ* nentries : " << data->GetEntries() << endl;
  cout << endl;

  //------Declaration of variables--------

  Int_t Nentries=data->GetEntries();
  TObjArray* array=0;
  NcTrack* tr=0;
  Nc3Vector Pos;
  NcTrack* Primtr=0;

  vector<Float_t> OAHMB;
  OAHMB.resize(NrecoHMB);
  vector<Float_t> OAB;
  OAB.resize(NrecoB);
  vector<Float_t> OARL;
  OARL.resize(NrecoRL);
  vector<Float_t> thetaHMB;
  thetaHMB.resize(NrecoHMB);
  vector<Float_t> thetaB;
  thetaB.resize(NrecoB);
  vector<Float_t> thetaRL;
  thetaRL.resize(NrecoRL);
  vector<Float_t> QualHMB;
  QualHMB.resize(NrecoHMB);
  vector<Float_t> QualB;
  QualB.resize(NrecoB);
  vector<Float_t> QualRL;
  QualRL.resize(NrecoRL);
  vector<Int_t> filter;
  filter.resize(7);


  Int_t c=0;
  Int_t d=0;
  Float_t theta=0.;
  Float_t beta=0.;
  Float_t RLogL=0.;
  NcSignal* fits=0;
  Double_t* DifZ=new Double_t;
  Double_t* Dist=new Double_t;
  Float_t condition=0.;
  NcDevice* filt=0;
  NcSignal* sx=0;
  Int_t nfilt=0;
  TString name;
  Float_t weight=0.;
  vector<Float_t> Prim;
  Prim.resize(3);
  //Float_t phi=0.;
  Int_t oldEvtN=0;
  Int_t EvtN=0;

 

  //Create the file, the Tree and a few branches
  TFile *file=new TFile(OUTfilename,"recreate");
  TTree *tree=new TTree("Data","Info about several tracks");
  
  //Reco info
  for(Int_t i=0;i<NrecoHMB;i++){
    tree->Branch(recoHMB[i]+"_theta",&thetaHMB[i],"theta/F");
    tree->Branch(recoHMB[i]+"_Qual",&QualHMB[i],"Qual/F");
    if(type!=2)tree->Branch(recoHMB[i]+"_OA",&OAHMB[i],"OA/F");
  }
  for(Int_t i=0;i<NrecoB;i++){
    tree->Branch(recoB[i]+"_theta",&thetaB[i],"theta/F");
    tree->Branch(recoB[i]+"_Qual",&QualB[i],"Qual/F");
    if(type!=2)tree->Branch(recoB[i]+"_OA",&OAB[i],"OA/F");
  }
  for(Int_t i=0;i<NrecoRL;i++){
    tree->Branch(recoRL[i]+"_theta",&thetaRL[i],"theta/F");
    tree->Branch(recoRL[i]+"_Qual",&QualRL[i],"Qual/F");
    if(type!=2)tree->Branch(recoRL[i]+"_OA",&OARL[i],"OA/F");
  }
  


  //Event info
  tree->Branch("Primary_energy",&Prim[0],"Primary_energy/F");
  tree->Branch("oneweight",&weight,"oneweight/F");
  tree->Branch("EvtN",&EvtN,"EvtN/I");
  tree->Branch("Primary_theta",&Prim[1],"Primary_theta/F");
  tree->Branch("Primary_phi",&Prim[2],"Primary_phi/F");
  
  //Filter info
  tree->Branch("ICOnlineL2",&filter[0],"ICOnlineL2/I");
  tree->Branch("SDST_MuonFilter",&filter[1],"SDST_MuonFilter/I");
  tree->Branch("MuonFilter",&filter[2],"MuonFilter/I");
  tree->Branch("LowUp",&filter[3],"LowUp/I");
  tree->Branch("EHE",&filter[4],"EHE/I");
  tree->Branch("Cascade",&filter[5],"Cascade/I");
  tree->Branch("DeepCore",&filter[6],"DeepCore/I");

  //Loop over all events
  for(Int_t ien=0;ien<Nentries;ien++)
    {
   
      //-----Get data------
      data->GetEntry(ien);
      if(!evt)continue;

      //-----Follow the evolution------
      if(c*100./Nentries>d)
	{
	  cout<<d<<"%"<<endl;
	  d++;
	}
      c++;


      //Check for non redundancy	 
      oldEvtN=EvtN;
      EvtN=evt->GetEventNumber();
      
      if(EvtN==oldEvtN && type==2)continue;
      
      //------Primary------
      //Store the primary info into Primary track                                                                                                            
    
      if(type!=2){
	Primtr=(NcTrack*)evt->GetIdTrack(-1);
	if(!Primtr){
	  cout<<"No Primary in evt "<<ien<<" !!"<<endl;
	  continue;
	}
	Pos=Primtr->Get3Vector();
	theta=Pos.GetX(2,"sph","deg");
	Prim[1]=theta;
	Prim[2]=Pos.GetX(3,"sph","deg");
	if(type==0){
	  if(theta>85.0 || theta<5.0) continue;
	}
	
	Prim[0]=Primtr->GetEnergy();
	weight=evt->GetWeight();

	
	
      }                                                                                                                  
	
     

      //----First type of recos: HMBeta recos-----  

      for(Int_t i=0;i<NrecoHMB;i++){

	OAHMB[i]=-1.;
	thetaHMB[i]=-1.;
	QualHMB[i]=0.;
	
	if(i==0 || i==1)array=evt->GetTracks(recoHMB[i],1);
	else array=evt->GetTracks(recoHMB[i],0);
	if(array->GetEntries()!=0)	  
	  {
	    tr=(NcTrack*)array->At(0);
	    if(tr->HasVector()==1){
	      Pos=tr->Get3Vector();
	      thetaHMB[i]=Pos.GetX(2,"sph","deg");
		
	      if(type!=2) OAHMB[i]=Primtr->GetOpeningAngle(Pos,"deg");
	    
	      beta=0.;
	      if(i==0 || i==1)beta=Beta(evt,recoHMB[i],1, DifZ, Dist);
	      else beta=Beta(evt,recoHMB[i],0, DifZ, Dist);
	    
	      if(beta<=1. && beta>=0.) QualHMB[i]=beta;
	      if(beta>1. && beta<2.) QualHMB[i]=2.-beta;
	      if(beta>2.)QualHMB[i]=0.;

	   
	    
	    }
	    
	  }

	
      }
      


      //----Second type: Beta recos-----  

    
      for(Int_t i=0;i<NrecoB;i++){

	OAB[i]=-1.;
	thetaB[i]=-1.;
	QualB[i]=0.;

	if(i==0)array=evt->GetTracks(recoB[i],1);
	else array=evt->GetTracks(recoB[i],0);
	
	if(array->GetEntries()!=0)	  
	  {
	   
	    tr=(NcTrack*)array->At(0);
	    if(tr->HasVector()==1){
	      Pos=tr->Get3Vector();
	      thetaB[i]=Pos.GetX(2,"sph","deg");

	      if(type!=2) OAB[i]=Primtr->GetOpeningAngle(Pos,"deg");
	    
	      beta=-1.;
	      fits=(NcSignal*)tr->GetFitDetails();
	      if (fits)
		{
		  beta=fits->GetSignal("Beta");
		  if(i==2)beta=fits->GetSignal("Speed")/0.299792;
		}
	      if(beta>=0 && beta<=1) QualB[i]=beta;
	      if(beta>1. && beta<2.) QualB[i]=2.-beta;
	      if(beta>2.)QualB[i]=0.;
	    
          
	    }
	    
	  }


	
      }
      
      //----Third type: RLogL reco-----  
     
      for(Int_t i=0;i<NrecoRL;i++){

	OARL[i]=-1.;
	thetaRL[i]=-1.;
	QualRL[i]=0.;

	if(i==2)array=evt->GetTracks(recoRL[i],1);
	else array=evt->GetTracks(recoRL[i],0);
	
	if(array->GetEntries()!=0)	  
	  {
	   
	    tr=(NcTrack*)array->At(0);
	   
	    if(tr->HasVector()==1){
	      Pos=tr->Get3Vector();
	      thetaRL[i]=Pos.GetX(2,"sph","deg");
		
	      if(type!=2) OARL[i]=Primtr->GetOpeningAngle(Pos,"deg");
	    
	      RLogL=-1.;
	      fits=(NcSignal*)tr->GetFitDetails();
	      if (fits)
		{
		  RLogL=fits->GetSignal("RLogL");
		}
	      if(RLogL>=6.) QualRL[i]=6./RLogL;
	      if(RLogL>0. && RLogL<6.) QualRL[i]=1.;
	    
     
	    }
	      
	  }

      }
      

      //Filter info
      for(Int_t i=0;i<7;i++)filter[i]=0;
      
      filt=evt->GetDevice("Filter");
      nfilt=0;
      if (filt) nfilt=filt->GetNhits();
      
      for (Int_t ifilt=1; ifilt<=nfilt; ifilt++)
	{
	  sx=filt->GetHit(ifilt);
	  if (!sx) continue;

	  name=sx->GetName();
	  condition=sx->GetSignal("condition");
	  
	  if (condition)
	    {
	      if (name.Contains("ICOnlineL2")) filter[0]=1;
	      
	      if (name.Contains("MuonFilter"))
		{
		  if (name.Contains("SDST")) filter[1]=1;
		  else filter[2]=1;
		}
	      
	      if (name.Contains("LowUp")) filter[3]=1;
	      
	      if (name.Contains("EHE")) filter[4]=1;
	      
	      if (name.Contains("Cascade")) filter[5]=1;
	      
	      if (name.Contains("DeepCore")) filter[6]=1;
	      
	    }
	}
      

      
      tree->Fill();
      
    }

 
  
  tree->Write();
  
  file->Close();
  
 
  delete[] recoRL;
  delete[] recoHMB;
  delete[] recoB;
  delete DifZ;
  delete Dist;
  delete evt;
  delete data;
 
  delete file;
 

  cout << "End" << endl ;

  return 0;
}