Interface to FASTJET. (Rafael Diaz)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Oct 2007 10:39:58 +0000 (10:39 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Oct 2007 10:39:58 +0000 (10:39 +0000)
JETAN/AliFastJetFinder.cxx [new file with mode: 0644]
JETAN/AliFastJetFinder.h [new file with mode: 0644]
JETAN/AliFastJetHeader.cxx [new file with mode: 0644]
JETAN/AliFastJetHeader.h [new file with mode: 0644]
JETAN/JETANLinkDef.h
JETAN/libJETAN.pkg

diff --git a/JETAN/AliFastJetFinder.cxx b/JETAN/AliFastJetFinder.cxx
new file mode 100644 (file)
index 0000000..f2c40c4
--- /dev/null
@@ -0,0 +1,523 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+//---------------------------------------------------------------------
+// FastJet finder
+// interface to FastJet algorithm
+// Author: Rafael.Diaz.Valdes@cern.ch
+// kt using NlnN 
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TLorentzVector.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TArrayF.h>
+#include <TRandom.h>
+#include <TClonesArray.h>
+
+#include "AliFastJetFinder.h"
+#include "AliFastJetHeader.h"
+#include "AliJetReaderHeader.h"
+#include "AliJetReader.h"
+#include "AliJet.h"
+//for FastJet finder
+#include "FjPseudoJet.hh"
+#include "FjClusterSequence.hh"
+
+
+ClassImp(AliFastJetFinder);
+
+
+////////////////////////////////////////////////////////////////////////
+
+AliFastJetFinder::AliFastJetFinder():
+    AliJetFinder(),
+    fHeader(0x0),
+    fLego(0x0),
+    fLegoSignal(0x0)
+{
+  // Constructor
+}
+
+////////////////////////////////////////////////////////////////////////
+
+AliFastJetFinder::~AliFastJetFinder()
+
+{
+  // destructor
+}
+
+////////////////////////////////////////////////////////////////////////
+
+
+void AliFastJetFinder::FindJets()
+
+{
+  //create cells and jets array
+  // 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 particles 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];
+  Int_t*   injet = new Int_t[nIn];
+
+
+  //total energy in array
+  Float_t  EtTotal = 0.0;
+  Float_t  meanptCell = 0.0;
+  Float_t  sqptCell = 0.0;
+
+
+  // load input vectors in fLego
+  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]);
+      if(fReader->GetSignalFlag(i) == 1)
+        fLegoSignal->Fill(etaT[i], phiT[i], ptT[i]);
+      EtTotal= EtTotal+ptT[i];
+    }
+  }
+  fJets->SetNinput(nIn);
+
+  // add soft background fixed
+  Int_t nsoft = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
+  Float_t* ptRndm =  new Float_t[nsoft];
+  if(fHeader->AddSoftBackg()){
+    gRandom->RndmArray(nsoft,ptRndm);
+    for(Int_t isoft = 0; isoft < nsoft; isoft++){
+        Float_t ptsoft  = 0.005*ptRndm[isoft];
+        EtTotal= EtTotal+ptsoft;
+    }
+  }
+
+  if(EtTotal == 0){
+     delete enT;
+     delete ptT;
+     delete etaT;
+     delete phiT;
+     return;
+  }
+
+  //cell array
+  Float_t* etCell = new Float_t[90000];   //! Cell Energy // check enough space! *to be done*
+  Float_t* etaCell = new Float_t[90000];  //! Cell eta
+  Float_t* phiCell = new Float_t[90000];  //! Cell phi
+  Int_t*   jetflagCell = new Int_t[90000]; //! Cell flag for jets
+  Float_t* etsigCell = new Float_t[90000];   // signal in this cell
+
+  //jets array
+  Float_t* etaJet = new Float_t[200];
+  Float_t* phiJet = new Float_t[200];
+  Float_t* etJet  = new Float_t[200];
+  Float_t* etsigJet  = new Float_t[200]; //signal energy
+  Float_t* etallJet = new Float_t[200];  //total energy in jet area
+  Int_t*   ncellsJet = new Int_t[200];
+  memset(etaJet,0,sizeof(Float_t)*200);
+  memset(phiJet,0,sizeof(Float_t)*200);
+  memset(etJet,0,sizeof(Float_t)*200);
+  memset(etsigJet,0,sizeof(Float_t)*200);
+  memset(etallJet,0,sizeof(Float_t)*200);
+  memset(ncellsJet,0,sizeof(Int_t)*200);
+
+
+
+  // load cells arrays
+  Int_t nCell = 0;
+  TAxis* xaxis = fLego->GetXaxis();
+  TAxis* yaxis = fLego->GetYaxis();
+  Float_t e = 0.0; Float_t esig = 0.0;
+
+  for (Int_t ib = 1; ib <= fHeader->GetLegoNbinEta(); ib++) {
+      for (Int_t jb = 1; jb <= fHeader->GetLegoNbinPhi(); jb++) {
+              e = fLego->GetBinContent(ib,jb);
+              if (e < 0.0) continue; // don't include this cells
+              Float_t eta  = xaxis->GetBinCenter(ib);
+              Float_t phi  = yaxis->GetBinCenter(jb);
+          if(fHeader->AddSoftBackg())
+             etCell[nCell]  = e + 0.005*ptRndm[nCell];
+          else
+             etCell[nCell]  = e;
+          sqptCell = sqptCell + etCell[nCell]*etCell[nCell]; // xi^2 ////////
+              etaCell[nCell] = eta;
+              phiCell[nCell] = phi;
+          jetflagCell[nCell] = -1; //default
+          esig = fLegoSignal->GetBinContent(ib,jb);
+          if(esig > 0.0)
+             etsigCell[nCell] = esig;
+          else
+             etsigCell[nCell] = 0.0;
+              nCell++;
+      }
+  }
+
+  meanptCell = EtTotal/(Float_t)nCell;
+  sqptCell = sqptCell/(Float_t)nCell;
+
+  Int_t nJets = 0;
+  //call to FastJet Algorithm
+  RunAlgorithm(nJets,etJet,etaJet,phiJet,etsigJet,etallJet,ncellsJet,
+               nCell,etCell,etaCell,phiCell,etsigCell,jetflagCell);
+
+
+  //subtract background
+  SubtractBackg(nCell,jetflagCell,etCell,
+                nJets,etJet,etallJet,ncellsJet,
+                meanptCell,sqptCell,EtTotal);
+
+
+  // add jets to list
+  Int_t* index = new Int_t[nJets];
+  Int_t nj = 0;
+  for(Int_t kj=0; kj<nJets; kj++){
+      if ((etaJet[kj] > (fHeader->GetJetEtaMax())) ||
+          (etaJet[kj] < (fHeader->GetJetEtaMin())) ||
+          (etJet[kj] < fHeader->GetMinJetEt())) continue; // acceptance eta range and etmin
+      Float_t px, py,pz,en; // convert to 4-vector
+      px = etJet[kj] * TMath::Cos(phiJet[kj]);
+      py = etJet[kj] * TMath::Sin(phiJet[kj]);
+      pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj])));
+      en = TMath::Sqrt(px * px + py * py + pz * pz);
+      fJets->AddJet(px, py, pz, en);
+      index[nj] = kj;
+      nj++;
+  }
+
+  //add signal percentage and total signal  in AliJets for analysis tool
+  Float_t* percentage  = new Float_t[nj];
+  Int_t* ncells      = new Int_t[nj];
+  Int_t* mult        = new Int_t[nj];
+
+  for(Int_t i = 0; i< nj; i++){
+     percentage[i] = etsigJet[index[i]]/etJet[index[i]];
+     ncells[i] = ncellsJet[index[i]];
+  }
+
+   //reorder injet flags
+  for(Int_t ipar = 0; ipar < nIn; ipar++){
+     Float_t injetflag =0;
+     Int_t iparCell = fLego->FindBin(etaT[ipar], phiT[ipar]);
+     injet[ipar] = jetflagCell[iparCell];
+     for(Int_t js = 0; js < nj; js++){
+        if(injet[ipar] == index[js]){
+          injet[ipar] = js;  // set the new jet id value
+          mult[js]++; // add multiplicity in jet js
+          injetflag = 1;
+          break;
+        }
+     }
+     if(injetflag == 0) injet[ipar] = -1; // set default value
+  }
+
+
+  fJets->SetNCells(ncells);
+  fJets->SetPtFromSignal(percentage);
+  fJets->SetMultiplicities(mult);
+  fJets->SetInJet(injet);
+  fJets->SetEtaIn(etaT);
+  fJets->SetPhiIn(phiT);
+  fJets->SetPtIn(ptT);
+  fJets->SetEtAvg(meanptCell);
+
+   //delete
+  delete enT;
+  delete ptT;
+  delete etaT;
+  delete phiT;
+  delete injet;
+  //cells
+  delete etCell;
+  delete etaCell;
+  delete phiCell;
+  delete jetflagCell;
+  delete etsigCell;
+  //jets
+  delete etaJet;
+  delete phiJet;
+  delete etJet;
+  delete etsigJet;
+  delete etallJet;
+  delete ncellsJet;
+
+  delete index;
+  delete percentage;
+  delete ncells;
+  delete mult;
+  delete ptRndm;
+
+}
+
+////////////////////////////////////////////////////////////////////////
+void AliFastJetFinder::RunAlgorithm(Int_t& nJets,Float_t* etJet,Float_t* etaJet,Float_t* phiJet,
+                                    Float_t* etsigJet, Float_t* etallJet, Int_t* ncellsJet,
+                                    Int_t& nCell,Float_t* etCell,Float_t* etaCell,Float_t* phiCell,
+                                    Float_t* etsigCell, Int_t* jetflagCell)
+{
+   //FastJet objects
+   vector<FjPseudoJet> input_cells; // create a vector
+   for (Int_t i = 0; i < nCell; i++){
+      if(etCell[i] == 0.0) continue; // not include cell empty
+      Double_t px, py,pz,en; // convert to 4-vector
+      px = etCell[i]*TMath::Cos(phiCell[i]);
+      py = etCell[i]*TMath::Sin(phiCell[i]);
+      pz = etCell[i]/TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaCell[i])));
+      en = TMath::Sqrt(px * px + py * py + pz * pz);
+      FjPseudoJet input_cell(px,py,pz,en); // create FjPseudoJet object
+      input_cell.set_user_index(i); //label the cell into Fastjet algortihm
+      //push FjPseudoJet of (px,py,pz,en) onto back of the input_cells
+      input_cells.push_back(input_cell);
+   }
+
+   //run the jet clustering with option R=1.0 and strategy= Best
+   Double_t Rparam = fHeader->GetRadius(); // default 1.0;
+   FjClusterSequence clust_seq(input_cells,Rparam);
+
+
+   //vector to get clusters
+   vector<FjPseudoJet> clusters;
+
+   ///////////////////////////////////////////////////////////////////////////
+   //extract the inclusive jets with pt> ptmin sorted by pt
+   //Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
+  // Double_t ptbgRc = EtBackgT*(Rparam*Rparam*TMath::Pi()/areaT);
+  // Double_t ptbgRcfluct = dEtTotal*Rparam*TMath::Sqrt(TMath::Pi()/areaT);
+  // Double_t ptmin = ptbgRc + ptbgRcfluct;
+   clusters = clust_seq.inclusive_jets(0);
+   //////////////////////////////////////////////////////////////////////////
+
+   /////////////////////////////////////////////////////////////////////////
+   //extract the exclusive jets with dcut = 25 GeV**2 and sort them in order
+   //of increasing pt
+   //Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
+   //Double_t ptbgRc = EtBackgT*(Rparam*Rparam*TMath::Pi()/areaT);
+   //Double_t ptbgRcfluct = dEtTotal*Rparam*TMath::Sqrt(TMath::Pi()/areaT);
+   //Double_t ptmin = ptbgRc + ptbgRcfluct;
+   //Double_t ktbackg = (fHeader->GetKfactor())*ptmin;
+   //Double_t dcut = ktbackg*ktbackg;
+   //clusters = sorted_by_pt(clust_seq.exclusive_jets(dcut));
+   //clusters = sorted_by_pt(clust_seq.exclusive_jets(5));
+   /////////////////////////////////////////////////////////////////////////
+
+
+   //cout << " ktClusters " << clusters.size() << endl;
+   nJets = (Int_t)clusters.size();
+   ////////////////////////////////////////////////////////////////////////
+   // add all clusters to jet arrays
+   for(UInt_t ij = 0; ij < clusters.size(); ij++){
+       //constituents
+       vector<FjPseudoJet> constituents = clust_seq.constituents(clusters[ij]);
+       //fill jet array info
+       ncellsJet[ij] = (Int_t)constituents.size();
+       phiJet[ij] = clusters[ij].phi();
+       Float_t angle = TMath::ATan(clusters[ij].perp()/clusters[ij].pz());
+       angle = ((angle < 0) ? angle + TMath::Pi() : angle);
+       etaJet[ij] = - TMath::Log(TMath::Tan(angle/2.0));
+       etJet[ij] = clusters[ij].perp();
+       //get constituents cells
+       for(UInt_t jc = 0; jc < constituents.size(); jc++){ // loop for all cells in ij cluster
+           Int_t jcell = constituents[jc].user_index();
+           jetflagCell[jcell] = ij; //flag this cell for jet
+           etsigJet[ij] = etsigJet[ij] + etsigCell[jcell]; // add signal of this cell
+           etallJet[ij] = etallJet[ij] + etCell[jcell];   // add total of this cell
+       }
+   }
+
+}
+////////////////////////////////////////////////////////////////////////
+void AliFastJetFinder::SubtractBackg(Int_t& nCell, Int_t* jetflagCell, Float_t* etCell,
+                                     Int_t& nJets, Float_t* etJet, Float_t* etallJet, Int_t* ncellsJet,
+                                     Float_t& meanptCell, Float_t& sqptCell, Float_t& etBackg)
+{
+   // simplest method: subtract background from external region to jets
+
+   //tmp array to flag jets
+   Int_t flagJet[200];
+   //tmp array to flag jet-cell status
+   Int_t tmpjetflagCell[90000];
+
+   Float_t etBackgOld = 0;
+   Float_t prec  = fHeader->GetPrecBg();
+   Float_t bgprec = 1;
+
+   while(bgprec > prec){
+        //clear tmpjetflagCell
+        memset(tmpjetflagCell,-1,sizeof(Int_t)*90000); // init with -1 (all cells are background)
+        //clear flagjet
+        memset(flagJet,0,sizeof(Int_t)*200); // init with 0 (no flag jets)
+        // select clusters > meantmpCell
+        for(Int_t i = 0; i < nJets; i++){
+            Float_t iptcell = etallJet[i]/(Float_t)ncellsJet[i];
+            if(iptcell < meanptCell) continue; // cluster not selected
+            // convert tmp cell background to jet cell
+            for(Int_t ic = 0; ic < nCell; ic++){ //loop for all cells
+                if(jetflagCell[ic] != i) continue; // other cells
+                tmpjetflagCell[ic] = i; // convert to a jet cell
+            }
+            //load total energy in cluster
+            etJet[i] = etallJet[i];
+            flagJet[i] = 1; // flag jet
+       }
+       //subtract background
+       for(Int_t j = 0; j < nCell; j++){ // loop for all cells
+           Int_t idxjet = tmpjetflagCell[j];
+           if(idxjet == -1) continue; // background cell
+           if(idxjet == -2) continue; // background protected cell
+           etJet[idxjet] = etJet[idxjet] - meanptCell;
+       }
+       // evaluate background fluctuations (rms value)
+       Float_t rmsptCell = TMath::Sqrt(sqptCell - meanptCell*meanptCell);
+       //fake jets
+       for(Int_t k = 0; k < nJets; k++){
+           if(flagJet[k] != 1) continue; // only flaged jets
+           //if(etJet[k] > fHeader->GetMinJetEt()) continue;  // jet ok!!
+           if(etJet[k] > rmsptCell*ncellsJet[k]) continue;  // jet ok!!
+           //clear tmpjetflag in cells
+           for(Int_t kc = 0; kc < nCell; kc++){ //loop for all cells
+               if(tmpjetflagCell[kc] != k) continue; // other cells
+               tmpjetflagCell[kc] = -1; // convert to background tmp cell
+           }
+           // clear all previous jet flags
+           etJet[k] = 0;
+           flagJet[k] = 0;
+       }
+       // recalculate background
+       sqptCell = 0;
+       etBackgOld = etBackg;
+       etBackg = 0;
+       Int_t nCellBackg = 0;
+       for(Int_t l = 0; l < nCell; l++){ // loop for all cells
+          if(tmpjetflagCell[l] != -1) continue; // cell included in some jet or protected
+          nCellBackg++;
+          etBackg = etBackg + etCell[l]; //add cell to background
+          //calc sqptCell
+          sqptCell = sqptCell + etCell[l]*etCell[l];
+       }
+       if(nCellBackg){
+          meanptCell = etBackg/(Float_t)nCellBackg; // new pt cell mean value
+          sqptCell = sqptCell/(Float_t)nCellBackg;
+       }else{
+          meanptCell = 0;
+          sqptCell = 0;
+       }
+       // evaluate presicion values
+       if(etBackg)
+         bgprec = (etBackgOld - etBackg)/etBackg;
+       else
+         bgprec = 0;
+   }
+
+   // set etJet 0 for all clusters not flaged in order to
+   for(Int_t m = 0; m < nJets; m++){
+       if(flagJet[m] == 1) continue; // flaged jets
+       etJet[m] = 0; //others clusters
+   }
+
+
+}
+
+////////////////////////////////////////////////////////////////////////
+void AliFastJetFinder::SubtractBackgArea(Int_t& nCell, Int_t* jetflagCell,
+                                     Float_t* etCell,Int_t&nJets, Float_t* etJet, Float_t* etallJet)
+{
+   // area method: subtract background from external region to jets
+   // using fixed area pi*Rc2
+
+   // n cells contained in a cone Rc
+   Double_t Rc = fHeader->GetRadius();
+   Float_t nCellRc = Rc*Rc*TMath::Pi()/(0.015*0.015); // change in future !!!!
+   //tmp array to flag fake jets
+   Int_t fakeflagJet[100];
+   memset(fakeflagJet,0,sizeof(Int_t)*100); // init with 0 (no fake jets)
+   Int_t njfake = nJets;
+   while(njfake > 0){
+       //calc background per cell
+       Int_t nCellBackg = 0;
+       Float_t EtBackg = 0.0;
+       for(Int_t i = 0; i < nCell; i++){ // loop for all cells
+           if(jetflagCell[i] != -1) continue; // cell included in some jet
+           nCellBackg++;
+           EtBackg = EtBackg + etCell[i]; //add cell to background
+       }
+       //subtract background energy per jet
+       for(Int_t l = 0; l < nJets; l++){
+           if(fakeflagJet[l] == 1) continue; // fake jet
+           etJet[l] = etallJet[l] - nCellRc*EtBackg/(Float_t)nCellBackg;
+       }
+       //fake jets analysis
+       njfake = 0;
+       for(Int_t k = 0; k < nJets; k++){
+           if(fakeflagJet[k] == 1) continue; // fake jet
+           if(etJet[k] < fHeader->GetMinJetEt()){  // a new fake jet
+              //clear jet flag in cells
+              for(Int_t kc = 0; kc < nCell; kc++){ //loop for all cells
+                  Int_t kidx = jetflagCell[kc];
+                  if(kidx != k) continue; // other cells
+                  jetflagCell[kc] = -1; // convert to background cell
+              }
+              fakeflagJet[k] = 1; // mark as a fake jet
+              njfake++; //count fakes in this loop
+           }
+       }
+   }
+}
+
+////////////////////////////////////////////////////////////////////////
+
+
+
+void AliFastJetFinder::Reset()
+{
+  fLego->Reset();
+  fLegoSignal->Reset();
+  fJets->ClearJets();
+
+}
+////////////////////////////////////////////////////////////////////////
+
+void AliFastJetFinder::WriteJHeaderToFile()
+{
+  fOut->cd();
+  fHeader->Write();
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliFastJetFinder::Init()
+{
+  // initializes some variables
+   // book lego
+  fLego = new
+    TH2F("legoH","eta-phi",
+        fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
+        fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(),
+        fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
+  fLegoSignal = new
+    TH2F("legoSignalH","eta-phi signal",
+        fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
+        fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(),
+        fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
+
+}
diff --git a/JETAN/AliFastJetFinder.h b/JETAN/AliFastJetFinder.h
new file mode 100644 (file)
index 0000000..4d92749
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef ALIFASTJETFINDER_H
+#define ALIFASTJETFINDER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+//---------------------------------------------------------------------
+// Fast Jet finder algorithm interface
+// manages the search for jets
+// Author: Rafael.Diaz.Valdes@cern.ch
+// kt algorithm using NlnN
+//---------------------------------------------------------------------
+
+#include "AliJetFinder.h"
+class AliFastJetHeader;
+
+class AliFastJetFinder : public AliJetFinder
+{
+ public:
+
+  AliFastJetFinder();
+  ~AliFastJetFinder();
+
+  // getters
+
+  // setters
+  void SetJetHeader(AliFastJetHeader* h) {fHeader= h;}
+  // others
+  void FindJets();
+  void RunAlgorithm(Int_t& nJets,Float_t* etJet,Float_t* etaJet,Float_t* phiJet,
+                    Float_t* etsigJet, Float_t* etallJet, Int_t* ncellsJet,
+                    Int_t& nCell,Float_t* etCell,Float_t* etaCell,Float_t* phiCell,
+                    Float_t* etsigCell, Int_t* jetflagCell);
+  void SubtractBackg(Int_t& nCell, Int_t* jetflagCell, Float_t* etCell,
+                     Int_t& nJets, Float_t* etJet, Float_t* etallJet, Int_t* ncellsJet,
+                     Float_t& meanptCell, Float_t& sqptCell, Float_t& etBackg);
+  void SubtractBackgArea(Int_t& nCell, Int_t* jetflagCell, Float_t* etCell,
+                     Int_t& nJets, Float_t* etJet, Float_t* etallJet);
+  void Reset();
+  void Init();
+  void WriteJHeaderToFile();
+
+ protected:
+
+  AliFastJetHeader* fHeader;         // pointer to jet header
+  TH2F*             fLego;           //! Lego Histo
+  TH2F*             fLegoSignal;    // ! Lego histogram for signal
+
+  ClassDef(AliFastJetFinder,1)
+};
+
+#endif
diff --git a/JETAN/AliFastJetHeader.cxx b/JETAN/AliFastJetHeader.cxx
new file mode 100644 (file)
index 0000000..ecfe7f8
--- /dev/null
@@ -0,0 +1,54 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//---------------------------------------------------------------------
+//  FastJet header class
+// Stores parameters of particle algoritm
+// Author: Rafael.Diaz.Valdes@cern.ch
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TMath.h>
+#include "AliFastJetHeader.h"
+
+ClassImp(AliFastJetHeader)
+
+////////////////////////////////////////////////////////////////////////
+
+AliFastJetHeader::AliFastJetHeader():
+    AliJetHeader("AliFastJetHeader"),
+    fLegoNbinEta(60),
+    fLegoNbinPhi(210),
+    fLegoEtaMin(-0.9),
+    fLegoEtaMax(0.9),
+    fLegoPhiMin(0.),
+    fLegoPhiMax(2. * TMath::Pi()),
+    fRadius(1.0),
+    fMinJetEt(10.0),
+    fSoftBackg(kTRUE),
+    fPrecBg(0.035) 
+{
+  // Constructor
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliFastJetHeader::PrintParameters() const
+{
+  // prints out parameters of jet algorithm
+
+  cout << " FastJet algorithm " << endl;
+
+}
diff --git a/JETAN/AliFastJetHeader.h b/JETAN/AliFastJetHeader.h
new file mode 100644 (file)
index 0000000..ce5689a
--- /dev/null
@@ -0,0 +1,66 @@
+#ifndef ALIFASTJETHEADER_H
+#define ALIFASTJETHEADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// FastJet algorithm
+//
+// Author: Rafael.Diaz.Valdes@cern.ch
+//---------------------------------------------------------------------
+#include "AliJetHeader.h"
+
+class AliFastJetHeader : public AliJetHeader
+{
+ public:
+  AliFastJetHeader();
+  virtual ~AliFastJetHeader() { }
+
+  // Getters
+  Int_t   GetLegoNbinEta()   const  {return fLegoNbinEta;}
+  Int_t   GetLegoNbinPhi()   const  {return fLegoNbinPhi;}
+  Float_t GetLegoEtaMin()    const  {return fLegoEtaMin;}
+  Float_t GetLegoEtaMax()    const  {return fLegoEtaMax;}
+  Float_t GetLegoPhiMin()    const  {return fLegoPhiMin;}
+  Float_t GetLegoPhiMax()    const  {return fLegoPhiMax;}
+  Float_t GetRadius()        const  {return fRadius;}
+  Float_t GetMinJetEt()      const  {return fMinJetEt;}
+  Bool_t  AddSoftBackg()     const  {return fSoftBackg;}
+  Float_t GetPrecBg()        const  {return fPrecBg;}
+  // Setters
+  void SetLegoNbinEta(Int_t f) {fLegoNbinEta=f;}
+  void SetLegoNbinPhi(Int_t f) {fLegoNbinPhi=f;}
+  void SetLegoEtaMin(Float_t f) {fLegoEtaMin=f;}
+  void SetLegoEtaMax(Float_t f) {fLegoEtaMax=f;}
+  void SetLegoPhiMin(Float_t f) {fLegoPhiMin=f;}
+  void SetLegoPhiMax(Float_t f) {fLegoPhiMax=f;}
+  void SetRadius(Float_t f)     {fRadius = f;}
+  void SetMinJetEt(Float_t f) {fMinJetEt=f;}
+  void SetSoftBackg(Bool_t f) {fSoftBackg=f;}
+  void SetPrecBg(Float_t f) {fPrecBg=f;}
+  // others
+  void PrintParameters() const;
+
+ protected:
+  // parameters for legos
+  Int_t   fLegoNbinEta;         //! number of cells in eta
+  Int_t   fLegoNbinPhi;         //! number of cells in phi
+  Float_t fLegoEtaMin;          //! minimum eta
+  Float_t fLegoEtaMax;          //! maximum eta
+  Float_t fLegoPhiMin;          //! minimun phi
+  Float_t fLegoPhiMax;          //! maximum phi
+  //algorithm parameters
+  Float_t fRadius;              // R parameter  (fastjet)
+  Float_t fMinJetEt;            // Min Et of jet (geneal)
+  Bool_t  fSoftBackg;           // add soft background
+  Float_t fPrecBg;              // max value of change for BG (in %)
+
+
+  ClassDef(AliFastJetHeader,1)
+};
+#endif
index 0731d4f3ab2eb60023d31f5ebee3dab55cc5d81b..f94acd3b7c8e7d9ea262135a4a3bbe0eb776b530 100644 (file)
@@ -39,4 +39,9 @@
 #pragma        link C++ class AliJetMCReader;
 #pragma        link C++ class AliJetMCReaderHeader;
 #pragma        link C++ class AliJetHistos+;
+#ifdef WITHFASTJET
+#pragma        link C++ class AliFastJetFinder+;
+#pragma        link C++ class AliFastJetHeader+;
+#endif
+
 #endif
index 98100120f2573da0fa70ed7d37cd589936154fc5..01e1e809896877c7bb686dcc852123856902ceb9 100644 (file)
@@ -24,7 +24,17 @@ SRCS =       AliJet.cxx AliJetHeader.cxx \
        AliJetMCReaderHeader.cxx \
        AliJetHistos.cxx
 
+ifneq ($(FASTJET),)
+       SRCS+= AliFastJetFinder.cxx AliFastJetHeader.cxx 
+endif
+
 
 HDRS:= $(SRCS:.cxx=.h) 
 DHDR= JETANLinkDef.h
 
+ifneq ($(FASTJET),)
+       EINCLUDE:= $(FASTJET)/include
+       CINTFLAGS += -DWITHFASTJET
+endif
+
+