JETAN/AliUA1JetFinder.cxx

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