[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 <TFile.h>
#include <TH2F.h>
#include "AliUA1JetFinder.h"
#include "AliUA1JetHeader.h"
#include "UA1Common.h"
#include "AliJetReaderHeader.h"
#include "AliJetReader.h"
#include "AliJet.h"


ClassImp(AliUA1JetFinder)

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

AliUA1JetFinder::AliUA1JetFinder():
  fHeader(0x0),
  fLego(0x0)
{
  // Default constructor
}

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

AliUA1JetFinder::~AliUA1JetFinder()

{
  // destructor
}

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

#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::FindJetsTPC()

{
  // 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* enT  = new Float_t[nIn];
  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);
    enT[i]  = lv->Energy();
    ptT[i]  = lv->Pt();
    etaT[i] = lv->Eta();
    phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
    if (fReader->GetCutFlag(i) == 1) 
      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->GetLegoNbinEta(); i++) {
      for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); 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->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
  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);

  // sort jets
  Int_t * idx  = new Int_t[UA1JETS.njet];
  TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);

  // download jet info.   
  for(Int_t i = 0; i < UA1JETS.njet; i++) {
    // reject events outside acceptable eta range
    if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
	|| ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
	  continue;
      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
  // find also percentage of pt from pythia
  Int_t* injet = new Int_t[nIn];
  Int_t* sflag = new Int_t[nIn];
  for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
  Int_t* mult  = new Int_t[fJets->GetNJets()];
  Int_t* ncell  = new Int_t[fJets->GetNJets()];
  Float_t* percentage  = new Float_t[fJets->GetNJets()];

  for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
      Float_t pt_sig = 0.0;
      mult[i] = 0;
      ncell[i] = UA1JETS.ncellj[i];
      for (Int_t j = 0; j < nIn; j++) {
	  Float_t deta = etaT[j] - fJets->GetEta(i);
	  Float_t dphi = phiT[j] - fJets->GetPhi(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] == 0) {
	      injet[j] = -(i+1);
	      if ((fReader->GetCutFlag(j)) == 1 &&
		  (etaT[j] < fHeader->GetLegoEtaMax()) &&
		  (etaT[j] > fHeader->GetLegoEtaMin())) {
		  injet[j] = i+1;
		  mult[i]++;
		  if (fReader->GetSignalFlag(j) == 1) {
		    pt_sig+=ptT[j];
		    sflag[j]=1;
		  }
	      }
	  }
      }
      percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
	((Double_t) fJets->GetPt(i));    
  }
  fJets->SetNCells(ncell);
  fJets->SetPtFromSignal(percentage);
  fJets->SetMultiplicities(mult);
  fJets->SetInJet(injet);
  fJets->SetEtaIn(etaT);
  fJets->SetPhiIn(phiT);
  fJets->SetPtIn(ptT);
  fJets->SetEtAvg(UA1JETS.etavg);
  delete ncell;
  delete enT;
  delete ptT;
  delete etaT;
  delete phiT;
  delete injet;
  delete idx;
  delete mult;
  delete percentage;
}

void AliUA1JetFinder::FindJets()
{ // Find jets with TPC or EMCAL or TPC + EMCAL information
  // initialize event, look for jets, download jet info

  // 1) transform input to pt,eta,phi plus lego
    AliJetUnitArray *fUnit = fReader->GetUnitArray();
    Int_t nIn = fUnit->GetUnitEntries();
    Int_t fOption = fReader->GetReaderHeader()->GetDetector();
    Int_t fDebug = fReader->GetReaderHeader()->GetDebug();
    
    if(fDebug>1){ 
	printf("Inside FindJets function ! \n");
	printf("fOption : %d \n",fOption);
    }
    
    if(fDebug>10){ 
	cout << "fUnit : " << fUnit << endl;
	printf("nIn = fUnit->GetUnitEntries() : %d \n",nIn);
    }
    
    if(fDebug>10){
	printf("     -----------------------------------------------------------------\n");
	printf("     --- All inputs in fUnitArray ---\n");
	for(Int_t i=0; i<nIn ; i++){
	    if(fUnit[i].GetUnitEnergy()!=0){
		printf("     -----------------------------------------------------------------\n");
		cout << "     |  ID : " << fUnit[i].GetUnitID() << "  |  Eta : " << fUnit[i].GetUnitEta() << "  |  Phi : " << fUnit[i].GetUnitPhi() << "  |  Energy : " << fUnit[i].GetUnitEnergy() << " |" <<endl;
	    }
	}
    }
    
    if (nIn == 0) return;
    
    Int_t nCandidateCut = 0;
    Int_t nCandidateCut2 = 0;
    Int_t nCandidate = 0;
    for (Int_t i = 0; i < nIn; i++){
	if(fUnit[i].GetUnitEnergy()>0. && fUnit[i].GetUnitCutFlag()==1) {
	    // Number of good tracks/digits which have passed the pt cut
	    nCandidateCut += 1;
	}
    if(fUnit[i].GetUnitEnergy()>0.) {
	// Number of good tracks/digits without pt cut
	nCandidate += 1;
    }
    }
    
    if(fDebug>=3){
    cout << "nCandidate : " << nCandidate << endl;
    cout << "nCandidateCut : " << nCandidateCut << endl;
    cout << "nCandidateCut2 : " << nCandidateCut2 << endl;
    }
    
  // local arrays for input No Cuts
  // Both pt < ptMin and pt > ptMin
    Float_t* enT       = new Float_t[nCandidate];
    Float_t* ptT       = new Float_t[nCandidate];
    Float_t* etaT      = new Float_t[nCandidate];
    Float_t* phiT      = new Float_t[nCandidate];
    Float_t* detaT     = new Float_t[nCandidate];
    Float_t* dphiT     = new Float_t[nCandidate];
    Float_t* cFlagT    = new Float_t[nCandidate];
    Float_t* cClusterT = new Float_t[nCandidate];
    Float_t* idT       = new Float_t[nCandidate];
    Int_t loop1 = 0;
    
    fJets->SetNinput(nCandidate);
    
    if(fDebug>3){
	cout << "nCandidate : " << nCandidate << endl;
	//    cout << "nMultCandidate : " << nMultCandidate << endl;
    }
    
    Float_t *etCell = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
    Float_t *etaCell = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
    Float_t *phiCell = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
    
    // Information extracted from fUnitArray
    for(Int_t i=0; i<nIn; i++) 
    {
	if(fUnit[i].GetUnitCutFlag()==1){ 
	    etCell[i] = fUnit[i].GetUnitEnergy();
	    if (etCell[i] > 0.0) etCell[i] -= fHeader->GetMinCellEt();
	    if (etCell[i] < 0.0) etCell[i] = 0.;
	    etaCell[i] = fUnit[i].GetUnitEta();
	    phiCell[i] = fUnit[i].GetUnitPhi();
	}
	else {
	    etCell[i]  = 0.;
	    etaCell[i] = fUnit[i].GetUnitEta();
	    phiCell[i] = fUnit[i].GetUnitPhi();
	}
	
	if(fUnit[i].GetUnitEnergy()>0.){
	    ptT[loop1]   = fUnit[i].GetUnitEnergy();
	    enT[loop1]   = fUnit[i].GetUnitEnergy();
	    etaT[loop1]  = fUnit[i].GetUnitEta();
	    phiT[loop1]  = fUnit[i].GetUnitPhi();
	    detaT[loop1] = fUnit[i].GetUnitDeta();
	    dphiT[loop1] = fUnit[i].GetUnitDphi();
	    cFlagT[loop1]= fUnit[i].GetUnitCutFlag();
	    idT[loop1]   = fUnit[i].GetUnitID();
	    loop1++;
	}
    }
    
    
    if(fDebug > 40) // For comparison
    {
	for(Int_t j=0; j<nIn; j++) {
	    if(etCell[j]>0){
		cout << "etCell[" << j << "] : " << etCell[j] << endl;
		cout << "etaCell[" << j << "] : " << etaCell[j] << endl;
		cout << "phiCell[" << j << "] : " << phiCell[j] << endl;
	    }
	}
    }
    
    // Run the algo. Parameters from header
    //  Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
    Int_t nTot      = nIn;
    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(nIn, nTot, etCell, etaCell, phiCell, 
		   minmove, maxmove, mode, precbg, ierror);
    
    // sort jets
    Int_t * idx  = new Int_t[UA1JETS.njet];
    TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
    
    if(fDebug > 20)
    {
	for(Int_t i = 0; i < UA1JETS.njet; i++) 
	{
	    cout << "Number of jets found, UA1JETS.njet : " << UA1JETS.njet << endl;
	    cout << "UA1JETS.etj : " << UA1JETS.etj << endl;
	    cout << "idx[" << i << "] : " << idx[i] << endl;
	    cout << "UA1JETS.etaj[1][" << idx[i] << "] : " << UA1JETS.etaj[1][idx[i]] << endl;
	    cout << "UA1JETS.phij[1][" << idx[i] << "] : " << UA1JETS.phij[1][idx[i]] << endl;
	    cout << "UA1JETS.etj[" << idx[i] << "] : " << UA1JETS.etj[idx[i]] << endl;
	}
    }
    
    // download jet info.   
    for(Int_t i = 0; i < UA1JETS.njet; i++) {
	// reject events outside acceptable eta range
	if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
	    || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
	{
	    continue;
	}
	
	Float_t px, py,pz,en,pT; // 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);
	pT = TMath::Sqrt(px * px + py * py);
	
	fJets->AddJet(px, py, pz, en);
	
    }

    // find multiplicities and relationship jet-particle
    // find also percentage of pt from pythia
    
    Int_t* injet = new Int_t[nCandidate];
    Int_t* sflag = new Int_t[nCandidate];
    for (Int_t i = 0; i < nCandidate; i++) {injet[i]= 0;sflag[i]=0;}
    Int_t* mult  = new Int_t[fJets->GetNJets()];
    Int_t* ncell  = new Int_t[fJets->GetNJets()];
    Float_t* percentage  = new Float_t[fJets->GetNJets()];
    
    // Instead of using etaT below, it would be interesting to use the previous fUnitArray object
    // With the particle ID, it is possible to to have access to its physical properties and one can,
    // for example, set if the corresponding particle is inside or outside the jet with the flag 
    // kOutJet/kInJet, other possibilities...
    
    for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
	Float_t pt_sig = 0.0;
	mult[i] = 0;
	ncell[i] = UA1JETS.ncellj[i];
	for (Int_t j = 0; j < nCandidate; j++) {
	    Float_t deta = etaT[j] - fJets->GetEta(i);
	    Float_t dphi = phiT[j] - fJets->GetPhi(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] == 0) {
		injet[j] = -(i+1);
		if(cFlagT[j] == 1 &&
		   (etaT[j] < fHeader->GetLegoEtaMax()) &&
		   (etaT[j] > fHeader->GetLegoEtaMin())) {
		    injet[j] = i+1;
		    mult[i]++;
		    pt_sig+=enT[j];
		    sflag[j]=1;
		}
	    }
	    if(fDebug>10){
		cout << "mult[" << i << "] : " << mult[i] << endl;
		cout << "ncell[" << i << "] : " << ncell[i] << endl;
	    }
	}
	percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
	    ((Double_t) fJets->GetPt(i));    
    }
    
    fJets->SetNCells(ncell);
    fJets->SetPtFromSignal(percentage);
    fJets->SetMultiplicities(mult);
    fJets->SetInJet(injet);
    fJets->SetEtaIn(etaT);
    if(fDebug>10){
	for(Int_t i=0; i<nCandidate ; i++){
	    cout << "phiT[" << i << "] : " << phiT[i] << endl;
	    cout << "etaT[" << i << "] : " << etaT[i] << endl;
	}
    }
    fJets->SetPhiIn(phiT);
    fJets->SetPtIn(enT);
    fJets->SetEtAvg(UA1JETS.etavg);
    delete etCell;
    delete etaCell;
    delete phiCell;
    delete ncell;
    delete cFlagT;
    delete cClusterT;
    delete enT;
    delete ptT;
    delete etaT;
    delete phiT;
    delete detaT;
    delete dphiT;
    delete injet;
    delete idx;
    delete mult;
    delete percentage;
    
}


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

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->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
    /((Float_t) fHeader->GetLegoNbinEta());
  dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
    /((Float_t) fHeader->GetLegoNbinPhi());

  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->GetLegoNbinEta(), fHeader->GetLegoEtaMin(), 
	 fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(), 
	 fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
}
Commit	Line	Data
99e5fe42	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
99e5fe42	15
	16	//---------------------------------------------------------------------
	17	// UA1 Jet finder
	18	// manages the search for jets
	19	// Author: jgcn@mda.cinvestav.mx
	20	// (code adapted from EMCAL directory)
	21	//---------------------------------------------------------------------
	22
	23	#include <Riostream.h>
	24	#include <TLorentzVector.h>
83a444b1	25	#include <TFile.h>
99e5fe42	26	#include <TH2F.h>
	27	#include "AliUA1JetFinder.h"
	28	#include "AliUA1JetHeader.h"
	29	#include "UA1Common.h"
8011d399	30	#include "AliJetReaderHeader.h"
99e5fe42	31	#include "AliJetReader.h"
	32	#include "AliJet.h"
	33
	34
	35	ClassImp(AliUA1JetFinder)
	36
	37	////////////////////////////////////////////////////////////////////////
	38
1b7d5d7e	39	AliUA1JetFinder::AliUA1JetFinder():
	40	fHeader(0x0),
	41	fLego(0x0)
99e5fe42	42	{
1b7d5d7e	43	// Default constructor
99e5fe42	44	}
	45
	46	////////////////////////////////////////////////////////////////////////
	47
	48	AliUA1JetFinder::~AliUA1JetFinder()
	49
	50	{
99e5fe42	51	// destructor
99e5fe42	52	}
	53
	54	////////////////////////////////////////////////////////////////////////
	55
	56	#ifndef WIN32
	57	# define ua1_jet_finder ua1_jet_finder_
	58	# define hf1 hf1_
	59	# define type_of_call
	60
	61	#else
	62	# define ua1_jet_finder UA1_JET_FINDER
	63	# define hf1 HF1
	64	# define type_of_call _stdcall
	65	#endif
	66
	67	extern "C" void type_of_call
	68	ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
	69	Float_t etc[60000], Float_t etac[60000],
	70	Float_t phic[60000],
	71	Float_t& min_move, Float_t& max_move, Int_t& mode,
	72	Float_t& prec_bg, Int_t& ierror);
	73
	74	extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
	75
	76	////////////////////////////////////////////////////////////////////////
	77
b45b0c92	78	void AliUA1JetFinder::FindJetsTPC()
99e5fe42	79
	80	{
	81	// initialize event, look for jets, download jet info
	82
99e5fe42	83	// initialization in 2 steps
	84	// 1) transform input to pt,eta,phi plus lego
	85	// 2) dump lego
	86
	87	// 1) transform input to pt,eta,phi plus lego
	88	TClonesArray *lvArray = fReader->GetMomentumArray();
	89	Int_t nIn = lvArray->GetEntries();
	90	if (nIn == 0) return;
	91
	92	// local arrays for input
83a444b1	93	Float_t* enT = new Float_t[nIn];
99e5fe42	94	Float_t* ptT = new Float_t[nIn];
	95	Float_t* etaT = new Float_t[nIn];
	96	Float_t* phiT = new Float_t[nIn];
	97
	98	// load input vectors
	99	for (Int_t i = 0; i < nIn; i++){
83a444b1	100	TLorentzVector lv = (TLorentzVector) lvArray->At(i);
	101	enT[i] = lv->Energy();
	102	ptT[i] = lv->Pt();
	103	etaT[i] = lv->Eta();
	104	phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
	105	if (fReader->GetCutFlag(i) == 1)
99e5fe42	106	fLego->Fill(etaT[i], phiT[i], ptT[i]);
	107	}
	108	fJets->SetNinput(nIn);
	109
	110	// 2) dump lego
	111	// check enough space! to be done
	112	Float_t etCell[60000]; //! Cell Energy
	113	Float_t etaCell[60000]; //! Cell eta
	114	Float_t phiCell[60000]; //! Cell phi
	115
	116	Int_t nCell = 0;
	117	TAxis* xaxis = fLego->GetXaxis();
	118	TAxis* yaxis = fLego->GetYaxis();
	119	Float_t e = 0.0;
83a444b1	120	for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
83a444b1	121	for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
99e5fe42	122	e = fLego->GetBinContent(i,j);
	123	if (e > 0.0) e -= fHeader->GetMinCellEt();
	124	if (e < 0.0) e = 0.;
	125	Float_t eta = xaxis->GetBinCenter(i);
	126	Float_t phi = yaxis->GetBinCenter(j);
	127	etCell[nCell] = e;
	128	etaCell[nCell] = eta;
	129	phiCell[nCell] = phi;
	130	nCell++;
	131	}
	132	}
	133
	134	// run the algo. Parameters from header
83a444b1	135	Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
99e5fe42	136	Float_t minmove = fHeader->GetMinMove();
	137	Float_t maxmove = fHeader->GetMaxMove();
	138	Int_t mode = fHeader->GetMode();
	139	Float_t precbg = fHeader->GetPrecBg();
	140	Int_t ierror;
	141	ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell,
	142	minmove, maxmove, mode, precbg, ierror);
	143
99e5fe42	144	// sort jets
	145	Int_t * idx = new Int_t[UA1JETS.njet];
	146	TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
8011d399	147
8011d399	148	// download jet info.
8011d399	149	for(Int_t i = 0; i < UA1JETS.njet; i++) {
8011d399	150	// reject events outside acceptable eta range
83a444b1	151	if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
	152	\|\| ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
	153	continue;
	154	Float_t px, py,pz,en; // convert to 4-vector
	155	px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
	156	py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
	157	pz = UA1JETS.etj[idx[i]] /
	158	TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
	159	en = TMath::Sqrt(px * px + py * py + pz * pz);
	160	fJets->AddJet(px, py, pz, en);
8011d399	161	}
99e5fe42	162
8011d399	163	// find multiplicities and relationship jet-particle
8011d399	164	// find also percentage of pt from pythia
99e5fe42	165	Int_t* injet = new Int_t[nIn];
83a444b1	166	Int_t* sflag = new Int_t[nIn];
83a444b1	167	for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
8011d399	168	Int_t* mult = new Int_t[fJets->GetNJets()];
83a444b1	169	Int_t* ncell = new Int_t[fJets->GetNJets()];
8011d399	170	Float_t* percentage = new Float_t[fJets->GetNJets()];
	171
	172	for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
83a444b1	173	Float_t pt_sig = 0.0;
	174	mult[i] = 0;
	175	ncell[i] = UA1JETS.ncellj[i];
	176	for (Int_t j = 0; j < nIn; j++) {
	177	Float_t deta = etaT[j] - fJets->GetEta(i);
	178	Float_t dphi = phiT[j] - fJets->GetPhi(i);
	179	if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
	180	if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
	181	Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
	182	if (dr < fHeader->GetRadius() && injet[j] == 0) {
	183	injet[j] = -(i+1);
	184	if ((fReader->GetCutFlag(j)) == 1 &&
	185	(etaT[j] < fHeader->GetLegoEtaMax()) &&
	186	(etaT[j] > fHeader->GetLegoEtaMin())) {
	187	injet[j] = i+1;
	188	mult[i]++;
	189	if (fReader->GetSignalFlag(j) == 1) {
	190	pt_sig+=ptT[j];
	191	sflag[j]=1;
	192	}
	193	}
	194	}
99e5fe42	195	}
83a444b1	196	percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
83a444b1	197	((Double_t) fJets->GetPt(i));
99e5fe42	198	}
83a444b1	199	fJets->SetNCells(ncell);
8011d399	200	fJets->SetPtFromSignal(percentage);
99e5fe42	201	fJets->SetMultiplicities(mult);
99e5fe42	202	fJets->SetInJet(injet);
83a444b1	203	fJets->SetEtaIn(etaT);
	204	fJets->SetPhiIn(phiT);
	205	fJets->SetPtIn(ptT);
	206	fJets->SetEtAvg(UA1JETS.etavg);
	207	delete ncell;
	208	delete enT;
	209	delete ptT;
	210	delete etaT;
	211	delete phiT;
	212	delete injet;
	213	delete idx;
	214	delete mult;
	215	delete percentage;
99e5fe42	216	}
99e5fe42	217
b45b0c92	218	void AliUA1JetFinder::FindJets()
	219	{ // Find jets with TPC or EMCAL or TPC + EMCAL information
	220	// initialize event, look for jets, download jet info
	221
	222	// 1) transform input to pt,eta,phi plus lego
	223	AliJetUnitArray *fUnit = fReader->GetUnitArray();
	224	Int_t nIn = fUnit->GetUnitEntries();
	225	Int_t fOption = fReader->GetReaderHeader()->GetDetector();
	226	Int_t fDebug = fReader->GetReaderHeader()->GetDebug();
	227
	228	if(fDebug>1){
	229	printf("Inside FindJets function ! \n");
	230	printf("fOption : %d \n",fOption);
	231	}
	232
	233	if(fDebug>10){
	234	cout << "fUnit : " << fUnit << endl;
	235	printf("nIn = fUnit->GetUnitEntries() : %d \n",nIn);
	236	}
	237
	238	if(fDebug>10){
	239	printf(" -----------------------------------------------------------------\n");
	240	printf(" --- All inputs in fUnitArray ---\n");
	241	for(Int_t i=0; i<nIn ; i++){
	242	if(fUnit[i].GetUnitEnergy()!=0){
	243	printf(" -----------------------------------------------------------------\n");
	244	cout << " \| ID : " << fUnit[i].GetUnitID() << " \| Eta : " << fUnit[i].GetUnitEta() << " \| Phi : " << fUnit[i].GetUnitPhi() << " \| Energy : " << fUnit[i].GetUnitEnergy() << " \|" <<endl;
	245	}
	246	}
	247	}
	248
	249	if (nIn == 0) return;
	250
	251	Int_t nCandidateCut = 0;
	252	Int_t nCandidateCut2 = 0;
	253	Int_t nCandidate = 0;
	254	for (Int_t i = 0; i < nIn; i++){
	255	if(fUnit[i].GetUnitEnergy()>0. && fUnit[i].GetUnitCutFlag()==1) {
	256	// Number of good tracks/digits which have passed the pt cut
	257	nCandidateCut += 1;
	258	}
	259	if(fUnit[i].GetUnitEnergy()>0.) {
	260	// Number of good tracks/digits without pt cut
	261	nCandidate += 1;
	262	}
	263	}
	264
	265	if(fDebug>=3){
	266	cout << "nCandidate : " << nCandidate << endl;
	267	cout << "nCandidateCut : " << nCandidateCut << endl;
	268	cout << "nCandidateCut2 : " << nCandidateCut2 << endl;
	269	}
	270
	271	// local arrays for input No Cuts
	272	// Both pt < ptMin and pt > ptMin
	273	Float_t* enT = new Float_t[nCandidate];
	274	Float_t* ptT = new Float_t[nCandidate];
	275	Float_t* etaT = new Float_t[nCandidate];
	276	Float_t* phiT = new Float_t[nCandidate];
	277	Float_t* detaT = new Float_t[nCandidate];
	278	Float_t* dphiT = new Float_t[nCandidate];
	279	Float_t* cFlagT = new Float_t[nCandidate];
	280	Float_t* cClusterT = new Float_t[nCandidate];
	281	Float_t* idT = new Float_t[nCandidate];
282	Int_t loop1 = 0;
283
284	fJets->SetNinput(nCandidate);
285
286	if(fDebug>3){
287	cout << "nCandidate : " << nCandidate << endl;
288	// cout << "nMultCandidate : " << nMultCandidate << endl;
289	}
290
291	Float_t *etCell = new Float_t[nIn]; //! Cell Energy - Extracted from UnitArray
292	Float_t *etaCell = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
293	Float_t *phiCell = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
294
295	// Information extracted from fUnitArray
296	for(Int_t i=0; i<nIn; i++)
297	{
298	if(fUnit[i].GetUnitCutFlag()==1){
299	etCell[i] = fUnit[i].GetUnitEnergy();
300	if (etCell[i] > 0.0) etCell[i] -= fHeader->GetMinCellEt();
301	if (etCell[i] < 0.0) etCell[i] = 0.;
302	etaCell[i] = fUnit[i].GetUnitEta();
303	phiCell[i] = fUnit[i].GetUnitPhi();
304	}
305	else {
306	etCell[i] = 0.;
307	etaCell[i] = fUnit[i].GetUnitEta();
308	phiCell[i] = fUnit[i].GetUnitPhi();
309	}
310
311	if(fUnit[i].GetUnitEnergy()>0.){
312	ptT[loop1] = fUnit[i].GetUnitEnergy();
313	enT[loop1] = fUnit[i].GetUnitEnergy();
314	etaT[loop1] = fUnit[i].GetUnitEta();
315	phiT[loop1] = fUnit[i].GetUnitPhi();
316	detaT[loop1] = fUnit[i].GetUnitDeta();
317	dphiT[loop1] = fUnit[i].GetUnitDphi();
318	cFlagT[loop1]= fUnit[i].GetUnitCutFlag();
319	idT[loop1] = fUnit[i].GetUnitID();
320	loop1++;
321	}
322	}
323
324
325	if(fDebug > 40) // For comparison
326	{
327	for(Int_t j=0; j<nIn; j++) {
328	if(etCell[j]>0){
329	cout << "etCell[" << j << "] : " << etCell[j] << endl;
330	cout << "etaCell[" << j << "] : " << etaCell[j] << endl;
331	cout << "phiCell[" << j << "] : " << phiCell[j] << endl;
332	}
333	}
334	}
335
336	// Run the algo. Parameters from header
337	// Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
338	Int_t nTot = nIn;
339	Float_t minmove = fHeader->GetMinMove();
340	Float_t maxmove = fHeader->GetMaxMove();
341	Int_t mode = fHeader->GetMode();
342	Float_t precbg = fHeader->GetPrecBg();
343	Int_t ierror;
344
345	ua1_jet_finder(nIn, nTot, etCell, etaCell, phiCell,
346	minmove, maxmove, mode, precbg, ierror);
347
348	// sort jets
349	Int_t * idx = new Int_t[UA1JETS.njet];
350	TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
351
352	if(fDebug > 20)
353	{
354	for(Int_t i = 0; i < UA1JETS.njet; i++)
355	{
356	cout << "Number of jets found, UA1JETS.njet : " << UA1JETS.njet << endl;
357	cout << "UA1JETS.etj : " << UA1JETS.etj << endl;
358	cout << "idx[" << i << "] : " << idx[i] << endl;
359	cout << "UA1JETS.etaj[1][" << idx[i] << "] : " << UA1JETS.etaj[1][idx[i]] << endl;
360	cout << "UA1JETS.phij[1][" << idx[i] << "] : " << UA1JETS.phij[1][idx[i]] << endl;
361	cout << "UA1JETS.etj[" << idx[i] << "] : " << UA1JETS.etj[idx[i]] << endl;
362	}
363	}
364
365	// download jet info.
366	for(Int_t i = 0; i < UA1JETS.njet; i++) {
367	// reject events outside acceptable eta range
368	if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
369	\|\| ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
370	{
371	continue;
372	}
373
374	Float_t px, py,pz,en,pT; // convert to 4-vector
375	px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
376	py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
377	pz = UA1JETS.etj[idx[i]] /
378	TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
379	en = TMath::Sqrt(px * px + py * py + pz * pz);
380	pT = TMath::Sqrt(px * px + py * py);
381
382	fJets->AddJet(px, py, pz, en);
383
384	}
385
386	// find multiplicities and relationship jet-particle
387	// find also percentage of pt from pythia
388
389	Int_t* injet = new Int_t[nCandidate];
390	Int_t* sflag = new Int_t[nCandidate];
391	for (Int_t i = 0; i < nCandidate; i++) {injet[i]= 0;sflag[i]=0;}
392	Int_t* mult = new Int_t[fJets->GetNJets()];
393	Int_t* ncell = new Int_t[fJets->GetNJets()];
394	Float_t* percentage = new Float_t[fJets->GetNJets()];
395
396	// Instead of using etaT below, it would be interesting to use the previous fUnitArray object
397	// With the particle ID, it is possible to to have access to its physical properties and one can,
398	// for example, set if the corresponding particle is inside or outside the jet with the flag
399	// kOutJet/kInJet, other possibilities...
400
401	for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
402	Float_t pt_sig = 0.0;
403	mult[i] = 0;
404	ncell[i] = UA1JETS.ncellj[i];
405	for (Int_t j = 0; j < nCandidate; j++) {
406	Float_t deta = etaT[j] - fJets->GetEta(i);
407	Float_t dphi = phiT[j] - fJets->GetPhi(i);
408	if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
409	if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
410	Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
411	if (dr < fHeader->GetRadius() && injet[j] == 0) {
412	injet[j] = -(i+1);
413	if(cFlagT[j] == 1 &&
414	(etaT[j] < fHeader->GetLegoEtaMax()) &&
415	(etaT[j] > fHeader->GetLegoEtaMin())) {
416	injet[j] = i+1;
417	mult[i]++;
418	pt_sig+=enT[j];
419	sflag[j]=1;
420	}
421	}
422	if(fDebug>10){
423	cout << "mult[" << i << "] : " << mult[i] << endl;
424	cout << "ncell[" << i << "] : " << ncell[i] << endl;
425	}
426	}
427	percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
428	((Double_t) fJets->GetPt(i));
429	}
430
431	fJets->SetNCells(ncell);
432	fJets->SetPtFromSignal(percentage);
433	fJets->SetMultiplicities(mult);
434	fJets->SetInJet(injet);
435	fJets->SetEtaIn(etaT);
436	if(fDebug>10){
437	for(Int_t i=0; i<nCandidate ; i++){
438	cout << "phiT[" << i << "] : " << phiT[i] << endl;
439	cout << "etaT[" << i << "] : " << etaT[i] << endl;
440	}
441	}
442	fJets->SetPhiIn(phiT);
443	fJets->SetPtIn(enT);
444	fJets->SetEtAvg(UA1JETS.etavg);
445	delete etCell;
446	delete etaCell;
447	delete phiCell;
448	delete ncell;
449	delete cFlagT;
450	delete cClusterT;
451	delete enT;
452	delete ptT;
453	delete etaT;
454	delete phiT;
455	delete detaT;
456	delete dphiT;
457	delete injet;
458	delete idx;
459	delete mult;
460	delete percentage;
461
462	}
463
464
99e5fe42	465	////////////////////////////////////////////////////////////////////////
	466
	467	void AliUA1JetFinder::Reset()
99e5fe42	468	{
	469	fLego->Reset();
	470	fJets->ClearJets();
	471	}
	472
	473	////////////////////////////////////////////////////////////////////////
	474
	475	void AliUA1JetFinder::WriteJHeaderToFile()
	476	{
	477	fOut->cd();
	478	fHeader->Write();
	479	}
	480
99e5fe42	481	////////////////////////////////////////////////////////////////////////
	482
	483	void AliUA1JetFinder::Init()
	484	{
	485	// initializes some variables
	486	Float_t dEta, dPhi;
83a444b1	487	dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
	488	/((Float_t) fHeader->GetLegoNbinEta());
	489	dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
	490	/((Float_t) fHeader->GetLegoNbinPhi());
99e5fe42	491
	492	UA1CELL.etaCellSize = dEta;
	493	UA1CELL.phiCellSize = dPhi;
	494
	495	// jet parameters
	496	UA1PARA.coneRad = fHeader->GetRadius();
	497	UA1PARA.etSeed = fHeader->GetEtSeed();
	498	UA1PARA.ejMin = fHeader->GetMinJetEt();
	499	UA1PARA.etMin = fHeader->GetMinCellEt();
	500
	501	// book lego
	502	fLego = new
	503	TH2F("legoH","eta-phi",
83a444b1	504	fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
	505	fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
	506	fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
99e5fe42	507	}