First commit of new jet reconstruction and analysis code to be used for the

[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinder.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
   
  
 
//---------------------------------------------------------------------
// UA1 Jet finder 
// manages the search for jets 
// Author: jgcn@mda.cinvestav.mx
// (code adapted from EMCAL directory)
//---------------------------------------------------------------------

#include <Riostream.h>
#include <TLorentzVector.h>
#include <TH2F.h>
#include "AliUA1JetFinder.h"
#include "AliUA1JetHeader.h"
#include "UA1Common.h"
#include "AliJetReader.h"
#include "AliJet.h"


ClassImp(AliUA1JetFinder)

////////////////////////////////////////////////////////////////////////

AliUA1JetFinder::AliUA1JetFinder()

{
  //
  // Constructor
  //
  fHeader = 0x0;
  fLego   = 0x0;
}

////////////////////////////////////////////////////////////////////////

AliUA1JetFinder::~AliUA1JetFinder()

{
  //
  // destructor
  //

  delete fLego;            
  
  // reset and delete header
}

////////////////////////////////////////////////////////////////////////

#ifndef WIN32
# define ua1_jet_finder ua1_jet_finder_
# define hf1 hf1_
# define type_of_call

#else
# define ua1_jet_finder UA1_JET_FINDER
# define hf1 HF1
# define type_of_call _stdcall
#endif

extern "C" void type_of_call
ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
               Float_t  etc[60000],  Float_t etac[60000],
               Float_t  phic[60000],
               Float_t& min_move, Float_t& max_move, Int_t& mode,
               Float_t& prec_bg,  Int_t& ierror);

extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);

////////////////////////////////////////////////////////////////////////

void AliUA1JetFinder::FindJets()

{
  // initialize event, look for jets, download jet info


  // initialization in 2 steps
  // 1) transform input to pt,eta,phi plus lego
  // 2) dump lego
  
  // 1) transform input to pt,eta,phi plus lego
  TClonesArray *lvArray = fReader->GetMomentumArray();
  Int_t nIn =  lvArray->GetEntries();
  if (nIn == 0) return;
    
  // local arrays for input
  Float_t* ptT  = new Float_t[nIn];
  Float_t* etaT = new Float_t[nIn];
  Float_t* phiT = new Float_t[nIn];

  // load input vectors
  for (Int_t i = 0; i < nIn; i++){
      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
      ptT[i]  = lv->Pt();
      etaT[i] = lv->Eta();
      phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
      fLego->Fill(etaT[i], phiT[i], ptT[i]);
  }
  fJets->SetNinput(nIn);
  
  // 2) dump lego
  // check enough space! *to be done*
  Float_t etCell[60000];   //! Cell Energy
  Float_t etaCell[60000];  //! Cell eta
  Float_t phiCell[60000];  //! Cell phi

  Int_t nCell = 0;
  TAxis* xaxis = fLego->GetXaxis();
  TAxis* yaxis = fLego->GetYaxis();
  Float_t e = 0.0;
  for (Int_t i = 1; i <= fHeader->GetNbinEta(); i++) {
      for (Int_t j = 1; j <= fHeader->GetNbinPhi(); j++) {
	  e = fLego->GetBinContent(i,j);
	  if (e > 0.0) e -= fHeader->GetMinCellEt();
	  if (e < 0.0) e = 0.;
	  Float_t eta  = xaxis->GetBinCenter(i);
	  Float_t phi  = yaxis->GetBinCenter(j);	    
	  etCell[nCell]  = e;
	  etaCell[nCell] = eta;
	  phiCell[nCell] = phi;
	  nCell++;
      }
  }

  // run the algo. Parameters from header
  Int_t nTot      = (fHeader->GetNbinEta())*(fHeader->GetNbinPhi());
  Float_t minmove = fHeader->GetMinMove();
  Float_t maxmove = fHeader->GetMaxMove();
  Int_t mode      = fHeader->GetMode();
  Float_t precbg  = fHeader->GetPrecBg();
  Int_t ierror;
  ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell, 
		 minmove, maxmove, mode, precbg, ierror);

  // download jet info. 
  Int_t* mult  = new Int_t[UA1JETS.njet];
  // sort jets
  Int_t * idx  = new Int_t[UA1JETS.njet];
  TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
  
  Int_t* injet = new Int_t[nIn];
  for (Int_t i = 0; i < nIn; i++) injet[i]= -1;

  for(Int_t i = 0; i < UA1JETS.njet; i++) {
      Float_t px, py,pz,en; // convert to 4-vector
      px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
      py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
      pz = UA1JETS.etj[idx[i]] /
	   TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
      en = TMath::Sqrt(px * px + py * py + pz * pz);
      
      fJets->AddJet(px, py, pz, en);
      // find multiplicities and relationship jet-particle
      mult[i] = 0;
      for (Int_t j = 0; j < nIn; j++) {
	  Float_t deta = etaT[j] - UA1JETS.etaj[1][idx[i]];
	  Float_t dphi = phiT[j] - UA1JETS.phij[1][idx[i]];
	  if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
	  if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
	  Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
	  if (dr < fHeader->GetRadius() && injet[j] == -1) {
	      injet[j] = i;
	      mult[i]++;
	  }
      }
  }
  fJets->SetMultiplicities(mult);
  fJets->SetInJet(injet);
}

////////////////////////////////////////////////////////////////////////

void AliUA1JetFinder::Reset()

{
  fLego->Reset();
  fJets->ClearJets();
}

////////////////////////////////////////////////////////////////////////

void AliUA1JetFinder::WriteJHeaderToFile()
{
  fOut->cd();
  fHeader->Write();
}


////////////////////////////////////////////////////////////////////////

void AliUA1JetFinder::Init()
{
  // initializes some variables
  Float_t dEta, dPhi;
  dEta = (fHeader->GetEtaMax()-fHeader->GetEtaMin())
    /((Float_t) fHeader->GetNbinEta());
  dPhi = (fHeader->GetPhiMax()-fHeader->GetPhiMin())
    /((Float_t) fHeader->GetNbinPhi());

  UA1CELL.etaCellSize = dEta;
  UA1CELL.phiCellSize = dPhi;

  // jet parameters
  UA1PARA.coneRad = fHeader->GetRadius();
  UA1PARA.etSeed  = fHeader->GetEtSeed();
  UA1PARA.ejMin   = fHeader->GetMinJetEt();
  UA1PARA.etMin   = fHeader->GetMinCellEt();

  // book lego
  fLego = new 
    TH2F("legoH","eta-phi",
	 fHeader->GetNbinEta(), fHeader->GetEtaMin(), fHeader->GetEtaMax(),
	 fHeader->GetNbinPhi(), fHeader->GetPhiMin(), fHeader->GetPhiMax());
}