[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()

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

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

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::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* 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::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
	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	**************************************************************************/
	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>
	25	#include <TFile.h>
	26	#include <TH2F.h>
	27	#include "AliUA1JetFinder.h"
	28	#include "AliUA1JetHeader.h"
	29	#include "UA1Common.h"
	30	#include "AliJetReaderHeader.h"
	31	#include "AliJetReader.h"
	32	#include "AliJet.h"
	33
	34
	35	ClassImp(AliUA1JetFinder)
	36
	37	////////////////////////////////////////////////////////////////////////
	38
	39	AliUA1JetFinder::AliUA1JetFinder()
	40
	41	{
	42	// Constructor
	43	fHeader = 0x0;
	44	fLego = 0x0;
	45	}
	46
	47	////////////////////////////////////////////////////////////////////////
	48
	49	AliUA1JetFinder::~AliUA1JetFinder()
	50
	51	{
	52	// destructor
	53	}
	54
	55	////////////////////////////////////////////////////////////////////////
	56
	57	#ifndef WIN32
	58	# define ua1_jet_finder ua1_jet_finder_
	59	# define hf1 hf1_
	60	# define type_of_call
	61
	62	#else
	63	# define ua1_jet_finder UA1_JET_FINDER
	64	# define hf1 HF1
	65	# define type_of_call _stdcall
	66	#endif
	67
	68	extern "C" void type_of_call
	69	ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
	70	Float_t etc[60000], Float_t etac[60000],
	71	Float_t phic[60000],
	72	Float_t& min_move, Float_t& max_move, Int_t& mode,
	73	Float_t& prec_bg, Int_t& ierror);
	74
	75	extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
	76
	77	////////////////////////////////////////////////////////////////////////
	78
	79	void AliUA1JetFinder::FindJets()
	80
	81	{
	82	// initialize event, look for jets, download jet info
	83
	84	// initialization in 2 steps
	85	// 1) transform input to pt,eta,phi plus lego
	86	// 2) dump lego
	87
	88	// 1) transform input to pt,eta,phi plus lego
	89	TClonesArray *lvArray = fReader->GetMomentumArray();
	90	Int_t nIn = lvArray->GetEntries();
	91	if (nIn == 0) return;
	92
	93	// local arrays for input
	94	Float_t* enT = new Float_t[nIn];
	95	Float_t* ptT = new Float_t[nIn];
	96	Float_t* etaT = new Float_t[nIn];
	97	Float_t* phiT = new Float_t[nIn];
	98
	99	// load input vectors
	100	for (Int_t i = 0; i < nIn; i++){
	101	TLorentzVector lv = (TLorentzVector) lvArray->At(i);
	102	enT[i] = lv->Energy();
	103	ptT[i] = lv->Pt();
	104	etaT[i] = lv->Eta();
	105	phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
	106	if (fReader->GetCutFlag(i) == 1)
	107	fLego->Fill(etaT[i], phiT[i], ptT[i]);
	108	}
	109	fJets->SetNinput(nIn);
	110
	111	// 2) dump lego
	112	// check enough space! to be done
	113	Float_t etCell[60000]; //! Cell Energy
	114	Float_t etaCell[60000]; //! Cell eta
	115	Float_t phiCell[60000]; //! Cell phi
	116
	117	Int_t nCell = 0;
	118	TAxis* xaxis = fLego->GetXaxis();
	119	TAxis* yaxis = fLego->GetYaxis();
	120	Float_t e = 0.0;
	121	for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
	122	for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
	123	e = fLego->GetBinContent(i,j);
	124	if (e > 0.0) e -= fHeader->GetMinCellEt();
	125	if (e < 0.0) e = 0.;
	126	Float_t eta = xaxis->GetBinCenter(i);
	127	Float_t phi = yaxis->GetBinCenter(j);
	128	etCell[nCell] = e;
	129	etaCell[nCell] = eta;
	130	phiCell[nCell] = phi;
	131	nCell++;
	132	}
	133	}
	134
	135	// run the algo. Parameters from header
	136	Int_t nTot = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
	137	Float_t minmove = fHeader->GetMinMove();
	138	Float_t maxmove = fHeader->GetMaxMove();
	139	Int_t mode = fHeader->GetMode();
	140	Float_t precbg = fHeader->GetPrecBg();
	141	Int_t ierror;
	142	ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell,
	143	minmove, maxmove, mode, precbg, ierror);
	144
	145	// sort jets
	146	Int_t * idx = new Int_t[UA1JETS.njet];
	147	TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
	148
	149	// download jet info.
	150	for(Int_t i = 0; i < UA1JETS.njet; i++) {
	151	// reject events outside acceptable eta range
	152	if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
	153	\|\| ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
	154	continue;
	155	Float_t px, py,pz,en; // convert to 4-vector
	156	px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
	157	py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
	158	pz = UA1JETS.etj[idx[i]] /
	159	TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
	160	en = TMath::Sqrt(px * px + py * py + pz * pz);
	161	fJets->AddJet(px, py, pz, en);
	162	}
	163
	164	// find multiplicities and relationship jet-particle
	165	// find also percentage of pt from pythia
	166	Int_t* injet = new Int_t[nIn];
	167	Int_t* sflag = new Int_t[nIn];
	168	for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
	169	Int_t* mult = new Int_t[fJets->GetNJets()];
	170	Int_t* ncell = new Int_t[fJets->GetNJets()];
	171	Float_t* percentage = new Float_t[fJets->GetNJets()];
	172
	173	for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
	174	Float_t pt_sig = 0.0;
	175	mult[i] = 0;
	176	ncell[i] = UA1JETS.ncellj[i];
	177	for (Int_t j = 0; j < nIn; j++) {
	178	Float_t deta = etaT[j] - fJets->GetEta(i);
	179	Float_t dphi = phiT[j] - fJets->GetPhi(i);
	180	if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
	181	if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
	182	Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
	183	if (dr < fHeader->GetRadius() && injet[j] == 0) {
	184	injet[j] = -(i+1);
	185	if ((fReader->GetCutFlag(j)) == 1 &&
	186	(etaT[j] < fHeader->GetLegoEtaMax()) &&
	187	(etaT[j] > fHeader->GetLegoEtaMin())) {
	188	injet[j] = i+1;
	189	mult[i]++;
	190	if (fReader->GetSignalFlag(j) == 1) {
	191	pt_sig+=ptT[j];
	192	sflag[j]=1;
	193	}
	194	}
	195	}
	196	}
	197	percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
	198	((Double_t) fJets->GetPt(i));
	199	}
	200	fJets->SetNCells(ncell);
	201	fJets->SetPtFromSignal(percentage);
	202	fJets->SetMultiplicities(mult);
	203	fJets->SetInJet(injet);
	204	fJets->SetEtaIn(etaT);
	205	fJets->SetPhiIn(phiT);
	206	fJets->SetPtIn(ptT);
	207	fJets->SetEtAvg(UA1JETS.etavg);
	208	delete ncell;
	209	delete enT;
	210	delete ptT;
	211	delete etaT;
	212	delete phiT;
	213	delete injet;
	214	delete idx;
	215	delete mult;
	216	delete percentage;
	217	}
	218
	219	////////////////////////////////////////////////////////////////////////
	220
	221	void AliUA1JetFinder::Reset()
	222	{
	223	fLego->Reset();
	224	fJets->ClearJets();
	225	}
	226
	227	////////////////////////////////////////////////////////////////////////
	228
	229	void AliUA1JetFinder::WriteJHeaderToFile()
	230	{
	231	fOut->cd();
	232	fHeader->Write();
	233	}
	234
	235	////////////////////////////////////////////////////////////////////////
	236
	237	void AliUA1JetFinder::Init()
	238	{
	239	// initializes some variables
	240	Float_t dEta, dPhi;
	241	dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
	242	/((Float_t) fHeader->GetLegoNbinEta());
	243	dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
	244	/((Float_t) fHeader->GetLegoNbinPhi());
	245
	246	UA1CELL.etaCellSize = dEta;
	247	UA1CELL.phiCellSize = dPhi;
	248
	249	// jet parameters
	250	UA1PARA.coneRad = fHeader->GetRadius();
	251	UA1PARA.etSeed = fHeader->GetEtSeed();
	252	UA1PARA.ejMin = fHeader->GetMinJetEt();
	253	UA1PARA.etMin = fHeader->GetMinCellEt();
	254
	255	// book lego
	256	fLego = new
	257	TH2F("legoH","eta-phi",
	258	fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
	259	fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
	260	fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
	261	}