Updates needed for current FASTJET (Rafael Diaz Valdes).
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 22 Sep 2008 11:23:05 +0000 (11:23 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 22 Sep 2008 11:23:05 +0000 (11:23 +0000)
JETAN/AliFastJetFinder.cxx
JETAN/AliFastJetFinder.h
JETAN/AliFastJetHeader.cxx
JETAN/AliFastJetHeader.h

index f2c40c4..6cd3e6e 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
  
+
 //---------------------------------------------------------------------
-// FastJet finder
-// interface to FastJet algorithm
+// FastJet v2.3.4 finder algorithm interface
+//
 // Author: Rafael.Diaz.Valdes@cern.ch
-// kt using NlnN 
+//  
 //---------------------------------------------------------------------
 
 #include <Riostream.h>
 #include "AliJetReaderHeader.h"
 #include "AliJetReader.h"
 #include "AliJet.h"
-//for FastJet finder
-#include "FjPseudoJet.hh"
-#include "FjClusterSequence.hh"
+
+
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/JetDefinition.hh"
+// get info on how fastjet was configured
+#include "fastjet/config.h"
+
+#ifdef ENABLE_PLUGIN_SISCONE
+#include "fastjet/SISConePlugin.hh"
+#endif
+
+#include<sstream>  // needed for internal io
+#include<vector> 
+#include <cmath> 
+
+using namespace std;
 
 
 ClassImp(AliFastJetFinder);
 
 
-////////////////////////////////////////////////////////////////////////
+//____________________________________________________________________________
 
 AliFastJetFinder::AliFastJetFinder():
-    AliJetFinder(),
-    fHeader(0x0),
-    fLego(0x0),
-    fLegoSignal(0x0)
+    AliJetFinder()
 {
   // Constructor
 }
 
-////////////////////////////////////////////////////////////////////////
+//____________________________________________________________________________
 
 AliFastJetFinder::~AliFastJetFinder()
-
 {
   // destructor
 }
 
-////////////////////////////////////////////////////////////////////////
-
-
+//______________________________________________________________________________
 void AliFastJetFinder::FindJets()
-
 {
-  //create cells and jets array
-  // 1) transform input to pt,eta,phi plus lego
+  
+  Bool_t debug = kFALSE;
+  
+  //pick up fastjet header
+  AliFastJetHeader *header = (AliFastJetHeader*)fHeader;
+
+  // check if we are reading AOD jets
+  TRefArray *refs = 0;
+  Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
+  if (fromAod) { refs = fReader->GetReferences(); }
+  
+  // RUN ALGORITHM  
+  // read input particles -----------------------------
+  vector<fastjet::PseudoJet> input_particles;
   TClonesArray *lvArray = fReader->GetMomentumArray();
+  if(lvArray == 0) { cout << "Could not get the momentum array" << endl; return; }
   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
+  if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }
+  Int_t nJets = 0; // n jets in this event
+  fJets->SetNinput(nIn) ; // number of input objects
+  Float_t px,py,pz,en;
+  // load input vectors
+  for(Int_t i = 0; i < nIn; i++){ // loop for all input particles
+      TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+      px = lv->Px();
+      py = lv->Py();
+      pz = lv->Pz();
+      en = lv->Energy();
+      
+      fastjet::PseudoJet input_part(px,py,pz,en); // create PseudoJet object
+      input_part.set_user_index(i); //label the particle into Fastjet algortihm
+      input_particles.push_back(input_part);  // back of the input_particles vector  
+  } // end loop 
+  
+  
+  // create an object that represents your choice of jet algorithm, and 
+  // the associated parameters
+  double Rparam = header->GetRparam();
+  fastjet::Strategy strategy = header->GetStrategy();
+  fastjet::RecombinationScheme recomb_scheme = header->GetRecombScheme();
+  fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); 
+  fastjet::JetDefinition jet_def(algorithm, Rparam, recomb_scheme, strategy);
+  
+  // create an object that specifies how we to define the area
+  fastjet::AreaDefinition area_def;
+  double ghost_etamax = header->GetGhostEtaMax(); 
+  double ghost_area   = header->GetGhostArea(); 
+  int    active_area_repeats = header->GetActiveAreaRepeats(); 
+  
+  // now create the object that holds info about ghosts
+  fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, ghost_area);
+  // and from that get an area definition
+  fastjet::AreaType area_type = header->GetAreaType();
+  area_def = fastjet::AreaDefinition(area_type,ghost_spec);
+  
+  // run the jet clustering with the above jet definition
+  fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def);
+  
+  
+  // save a comment in the header
+  TString comment = "Running FastJet algorithm with the following setup. ";
+  comment+= "Jet definition: ";
+  comment+= TString(jet_def.description());
+  comment+= ". Area definition: ";
+  comment+= TString(area_def.description());
+  comment+= ". Strategy adopted by FastJet: ";
+  comment+= TString(clust_seq.strategy_string());
+  header->SetComment(comment);
+  if(debug){
+    cout << "--------------------------------------------------------" << endl;
+    cout << comment << endl;
+    cout << "--------------------------------------------------------" << endl;
   }
-
-
-  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;
-
+  //header->PrintParameters();
+  
+  
+  // extract the inclusive jets with pt > ptmin, sorted by pt
+  double ptmin = header->GetPtMin(); 
+  vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
+  
+  //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
+  //subtract background // ===========================================
+  // set the rapididty , phi range within which to study the background 
+  double rap_max = header->GetRapMax(); 
+  double rap_min = header->GetRapMin();
+  double phi_max = header->GetPhiMax();
+  double phi_min = header->GetPhiMin();
+  fastjet::RangeDefinition range(rap_min, rap_max, phi_min, phi_max);
+  // subtract background
+  vector<fastjet::PseudoJet> sub_jets =  clust_seq.subtracted_jets(range,ptmin);  
+  
+  // print out
+  //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n";
+  //cout << "---------------------------------------\n";
+  //cout << endl;
+  //printf(" ijet   rap      phi        Pt         area  +-   err\n");
+   
+  // sort jets into increasing pt
+  vector<fastjet::PseudoJet> jets = sorted_by_pt(sub_jets);  
+  for (size_t j = 0; j < jets.size(); j++) { // loop for jets
+
+    double area     = clust_seq.area(jets[j]);
+    double area_error = clust_seq.area_error(jets[j]);
+
+    //printf("%5u %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, area_error);
+       
+       // go to write AOD  info
+    AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
+   //cout << "Printing jet " << endl;
+    if(debug) aodjet.Print("");
+   //cout << "Adding jet ... " ;
+    AddJet(aodjet);
+   //cout << "added \n" << endl;
+   
+  } // end loop for jets
+
+    
 }
 
-////////////////////////////////////////////////////////////////////////
-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::RunTest(const char* datafile)
 
-}
-////////////////////////////////////////////////////////////////////////
-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
+   // This simple test run the kt algorithm for an ascii file testdata.dat
+   // read input particles -----------------------------
+  vector<fastjet::PseudoJet> input_particles;
+  Float_t px,py,pz,en;
+  ifstream in;
+  Int_t nlines = 0;
+  // we assume a file basic.dat in the current directory
+  // this file has 3 columns of float data
+  in.open(datafile);
+  while (1) {
+      in >> px >> py >> pz >> en;
+      if (!in.good()) break;
+      //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz);
+      nlines++;
+      input_particles.push_back(fastjet::PseudoJet(px,py,pz,en)); 
    }
+   //printf(" found %d pointsn",nlines);
+   in.close();
+   //////////////////////////////////////////////////
+  // create an object that represents your choice of jet algorithm, and 
+  // the associated parameters
+  double Rparam = 1.0;
+  fastjet::Strategy strategy = fastjet::Best;
+  fastjet::RecombinationScheme recomb_scheme = fastjet::BIpt_scheme;
+  fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, recomb_scheme, strategy);
+  
+  
+  // create an object that specifies how we to define the area
+  fastjet::AreaDefinition area_def;
+  double ghost_etamax = 7.0;
+  double ghost_area    = 0.05;
+  int    active_area_repeats = 1;
+  
+
+  // now create the object that holds info about ghosts
+  fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, ghost_area);
+  // and from that get an area definition
+  area_def = fastjet::AreaDefinition(fastjet::active_area,ghost_spec);
+  
+
+  // run the jet clustering with the above jet definition
+  fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def);
+  
+  
+  // tell the user what was done
+  cout << "--------------------------------------------------------" << endl;
+  cout << "Jet definition was: " << jet_def.description() << endl;
+  cout << "Area definition was: " << area_def.description() << endl;
+  cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl;
+  cout << "--------------------------------------------------------" << endl;
+  
+  // extract the inclusive jets with pt > 5 GeV, sorted by pt
+  double ptmin = 5.0;
+  vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);
+  
+  cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
+  //subtract background // ===========================================
+  // set the rapididty range within which to study the background 
+  double rap_max = ghost_etamax - Rparam;
+  fastjet::RangeDefinition range(rap_max);
+  // subtract background
+  vector<fastjet::PseudoJet> sub_jets =  clust_seq.subtracted_jets(range,ptmin);  
+  
+  // print them out //================================================
+  cout << "Printing inclusive jets  after background subtraction \n";
+  cout << "------------------------------------------------------\n";
+  // sort jets into increasing pt
+  vector<fastjet::PseudoJet> jets = sorted_by_pt(sub_jets);  
+
+  printf(" ijet   rap      phi        Pt         area  +-   err\n");
+  for (size_t j = 0; j < jets.size(); j++) {
+
+    double area     = clust_seq.area(jets[j]);
+    double area_error = clust_seq.area_error(jets[j]);
+
+    printf("%5u %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",j,jets[j].rap(),
+          jets[j].phi(),jets[j].perp(), area, area_error);
+  }
+  cout << endl;
+  // ================================================================
 
-
-}
-
-////////////////////////////////////////////////////////////////////////
-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());
-
-}
+//____________________________________________________________________________
\ No newline at end of file
index 4d92749..55e6033 100644 (file)
@@ -5,15 +5,35 @@
  * See cxx source for full Copyright notice                               */
 
 
+
 //---------------------------------------------------------------------
-// Fast Jet finder algorithm interface
-// manages the search for jets
+// FastJet v2.3.4 finder algorithm interface
+//
 // Author: Rafael.Diaz.Valdes@cern.ch
-// kt algorithm using NlnN
+//  
 //---------------------------------------------------------------------
 
+//FastJet classes 
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/JetDefinition.hh"
+// get info on how fastjet was configured
+#include "fastjet/config.h"
+#ifdef ENABLE_PLUGIN_SISCONE
+#include "fastjet/SISConePlugin.hh"
+#endif
+
+
+#include<sstream>  // needed for internal io
+#include <vector> 
+#include <cmath> 
+
 #include "AliJetFinder.h"
-class AliFastJetHeader;
+#include "AliFastJetHeader.h"
+
+using namespace std;
+
 
 class AliFastJetFinder : public AliJetFinder
 {
@@ -22,32 +42,18 @@ class AliFastJetFinder : public AliJetFinder
   AliFastJetFinder();
   ~AliFastJetFinder();
 
-  // getters
-
-  // setters
-  void SetJetHeader(AliFastJetHeader* h) {fHeader= h;}
+  void FindJets(); 
   // 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 RunTest(const char* datafile); // a simple test
+  
   void WriteJHeaderToFile();
+  
+  protected:
+  AliFastJetFinder(const AliFastJetFinder& rfj);
+  AliFastJetFinder& operator = (const AliFastJetFinder& rsfj);
 
- protected:
-
-  AliFastJetHeader* fHeader;         // pointer to jet header
-  TH2F*             fLego;           //! Lego Histo
-  TH2F*             fLegoSignal;    // ! Lego histogram for signal
 
-  ClassDef(AliFastJetFinder,1)
+  ClassDef(AliFastJetFinder,2)
 };
 
 #endif
index ecfe7f8..92b9272 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+
 //---------------------------------------------------------------------
-//  FastJet header class
-// Stores parameters of particle algoritm
+// FastJet v2.3.4 finder algorithm interface
+// Finder Header Class 
 // Author: Rafael.Diaz.Valdes@cern.ch
 //---------------------------------------------------------------------
 
 #include <Riostream.h>
 #include <TMath.h>
+
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/JetDefinition.hh"
+
 #include "AliFastJetHeader.h"
 
 ClassImp(AliFastJetHeader)
@@ -29,18 +35,24 @@ 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) 
+    fRparam(1.0), 
+    fAlgorithm(fastjet::kt_algorithm),
+    fStrategy(fastjet::Best),
+    fRecombScheme(fastjet::BIpt_scheme),
+    fGhostEtaMax(2.0),
+    fGhostArea(0.05),
+    fActiveAreaRepeats(1),
+    fAreaType(fastjet::active_area), 
+    fPtMin(5.0),
+    fPhiMin(0),
+    fPhiMax(TMath::TwoPi())
 {
   // Constructor
+  
+  Double_t rapmax = fGhostEtaMax - fRparam;
+  Double_t rapmin = -fGhostEtaMax + fRparam;
+  SetRapRange(rapmin, rapmax);
+  
 }
 
 ////////////////////////////////////////////////////////////////////////
@@ -49,6 +61,27 @@ void AliFastJetHeader::PrintParameters() const
 {
   // prints out parameters of jet algorithm
 
-  cout << " FastJet algorithm " << endl;
+  cout << "FastJet algorithm  parameters:" << endl;
+  
+  cout << "-- Jet Definition --- " << endl;
+  cout << "R " << fRparam << endl;
+  cout << "Jet Algorithm " << fAlgorithm << endl; 
+  cout << "Strategy " << fStrategy << endl;  
+  cout << "Recombination Scheme " << fRecombScheme << endl; 
+  
+  cout << "-- Ghosted Area Spec parameters --- " << endl;
+  cout << "Ghost Eta Max " << fGhostEtaMax << endl;
+  cout << "Ghost Area " << fGhostArea << endl;
+  cout << "Active Area Repeats " << fActiveAreaRepeats << endl;
+  
+  cout << "-- Area Definition parameters --- " << endl;
+  cout << "Area Type " << fAreaType << endl; 
+  
+  cout << "-- Cluster Sequence Area parameters --- " << endl;
+  cout << "pt min " << fPtMin << endl; 
+  
+  cout << "-- Range Definition parameters --- " << endl;
+  cout << " bkg rapidity range from  " << fRapMin << " to " << fRapMax << endl;
+  cout << " bkg phi range from " << fPhiMin << " to " << fPhiMax << endl;
 
 }
index ce5689a..864c151 100644 (file)
@@ -5,11 +5,16 @@
  * See cxx source for full Copyright notice                               */
  
 //---------------------------------------------------------------------
-// FastJet algorithm
-//
+// FastJet v2.3.4 finder algorithm interface
+// Finder Header Class 
 // Author: Rafael.Diaz.Valdes@cern.ch
 //---------------------------------------------------------------------
  
+
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/JetDefinition.hh"
+
 #include "AliJetHeader.h"
 
  
@@ -21,46 +26,62 @@ class AliFastJetHeader : public AliJetHeader
   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;}
+  Double_t                     GetRparam()            const {return fRparam;}
+  fastjet::JetAlgorithm        GetAlgorithm()         const {return fAlgorithm;}
+  fastjet::Strategy            GetStrategy()          const {return fStrategy;}
+  fastjet::RecombinationScheme GetRecombScheme()      const {return fRecombScheme;}
+  Double_t                     GetGhostEtaMax()       const {return fGhostEtaMax;}
+  Double_t                     GetGhostArea()         const {return fGhostArea;}
+  Int_t                        GetActiveAreaRepeats() const {return fActiveAreaRepeats;}
+  fastjet::AreaType            GetAreaType()          const {return fAreaType;}
+  Double_t                     GetPtMin()             const {return fPtMin;}
+  Double_t                     GetRapMax()            const {return fRapMax;}
+  Double_t                     GetRapMin()            const {return fRapMin;}
+  Double_t                     GetPhiMax()            const {return fPhiMax;}
+  Double_t                     GetPhiMin()            const {return fPhiMin;}
+  
   // 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;}
+  void SetRparam(Double_t f)                           {fRparam = f;}
+  void SetAlgorithm(fastjet::JetAlgorithm f)           {fAlgorithm = f;}
+  void SetStrategy(fastjet::Strategy f)                {fStrategy = f;}
+  void SetRecombScheme(fastjet::RecombinationScheme f) {fRecombScheme = f;}
+  void SetGhostEtaMax(Double_t f)                      {fGhostEtaMax = f;}
+  void SetGhostArea(Double_t f)                        {fGhostArea = f;}
+  void SetActiveAreaRepeats(Int_t f)                   {fActiveAreaRepeats =f;}
+  void SetAreaType(fastjet::AreaType f)                {fAreaType = f;}
+  void SetRapRange(Double_t fmin, Double_t fmax)       {fRapMin = fmin; fRapMax = fmax;}
+  void SetPhiRange(Double_t fmin, Double_t fmax)       {fPhiMin = fmin; fPhiMax = fmax;}
+  
+  void SetComment(TString com) {fComment=com;}
+  
   // 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)
+  //fastjet::JetDefinition parameters
+  Double_t fRparam;
+  fastjet::JetAlgorithm fAlgorithm; //fastjet::kt_algorithm
+  fastjet::Strategy fStrategy;  //= fastjet::Best;
+  fastjet::RecombinationScheme fRecombScheme; // = fastjet::BIpt_scheme;
+  
+  //fastjet::GhostedAreaSpec parameters
+  Double_t fGhostEtaMax;
+  Double_t fGhostArea;
+  Int_t    fActiveAreaRepeats;
+  
+  //fastjet::AreaDefinition parameters
+  fastjet::AreaType fAreaType; 
+  
+  //fastjet::ClusterSequenceArea options parameters
+  Double_t fPtMin; //jets with pt > ptmin
+  
+  //fastjet::RangeDefinition parameters 
+  Double_t fRapMax, fRapMin; // rapidity range of background sub 
+  Double_t fPhiMax, fPhiMin; // phi range of background sub
+  
+  
+  ClassDef(AliFastJetHeader,2)
 };
  
 #endif