Interface to FASTJETs SISCone jetfinder added.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Apr 2009 13:23:55 +0000 (13:23 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 30 Apr 2009 13:23:55 +0000 (13:23 +0000)
JETAN/AliSISConeJetFinder.cxx [new file with mode: 0644]
JETAN/AliSISConeJetFinder.h [new file with mode: 0644]
JETAN/AliSISConeJetHeader.cxx [new file with mode: 0644]
JETAN/AliSISConeJetHeader.h [new file with mode: 0644]
JETAN/JETANLinkDef.h
JETAN/libJETAN.pkg

diff --git a/JETAN/AliSISConeJetFinder.cxx b/JETAN/AliSISConeJetFinder.cxx
new file mode 100644 (file)
index 0000000..307a2b7
--- /dev/null
@@ -0,0 +1,179 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//---------------------------------------------------------------------
+// SISCone (FastJet v2.3.4) finder algorithm interface
+//
+// Author: swensy.jangal@ires.in2p3.fr 
+//  
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TArrayF.h>
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TLorentzVector.h>
+#include <TRandom.h>
+
+#include "AliHeader.h"
+#include "AliJet.h"
+#include "AliJetKineReader.h"
+#include "AliJetReader.h"
+#include "AliJetReaderHeader.h"
+#include "AliSISConeJetFinder.h"
+#include "AliSISConeJetHeader.h"
+
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/PseudoJet.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(AliSISConeJetFinder)
+
+//____________________________________________________________________________
+
+AliSISConeJetFinder::AliSISConeJetFinder():
+AliJetFinder()
+{
+  // Constructor
+}
+
+//____________________________________________________________________________
+
+AliSISConeJetFinder::~AliSISConeJetFinder()
+{
+  // destructor
+}
+
+//______________________________________________________________________________
+void AliSISConeJetFinder::FindJets()
+{
+
+  Bool_t debug = kFALSE;
+  
+  // Pick up siscone header
+  AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader;
+
+  // Check if we are reading AOD jets
+  TRefArray *refs = 0;
+  Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
+  if (fromAod) { refs = fReader->GetReferences(); }
+  
+
+//******************************** SISCONE PLUGIN CONFIGURATION
+// Here we look for SISCone parameters in the header and we define our plugin.  
+
+  Double_t coneRadius                = header->GetConeRadius();                 // cone radius
+  Double_t overlapThreshold          = header->GetOverlapThreshold();           // overlap parameter
+  Int_t    nPassMax                  = header->GetNPassMax();                   // maximum number of passes
+  Double_t ptProtoJetMin             = header->GetPtProtojetMin();              // pT min of protojets
+  Double_t caching                   = header->GetCaching();                    // do we record found cones for this set of data?
+
+  if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
+
+  fastjet::JetDefinition::Plugin * plugin;
+  plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
+
+//********************************  READING OF INPUT PARTICLES
+// Here we look for px, py pz and energy of each particle that we gather in a PseudoJet object, and we put all these PseudoJet in a vector of PseudoJets : input_particles. 
+
+  TClonesArray *lvArray = fReader->GetMomentumArray();
+  Int_t nIn = lvArray->GetEntries();
+
+  // We check if lvArray is ok
+  if(lvArray == 0)
+  {
+    cout << "Could not get the momentum array" << endl;
+    return;
+  }
+
+  if(nIn == 0)// nIn = Number of particles in the event
+  {
+    if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
+    return;
+  }
+
+  Int_t nJets = 0;        // Number of jets in this event
+  fJets->SetNinput(nIn) ; // fJets = AliJet number of input objects
+  Float_t px,py,pz,en;
+  vector<fastjet::PseudoJet> inputParticles; 
+
+  // Load input vectors
+  for(Int_t i = 0; i < nIn; i++)
+  { 
+    TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
+    px = lv->Px();
+    py = lv->Py();
+    pz = lv->Pz();
+    en = lv->Energy();
+      
+    fastjet::PseudoJet inputPart(px,py,pz,en); 
+    inputPart.set_user_index(i);
+    inputParticles.push_back(inputPart); 
+
+  }   
+  
+//******************************** JETS FINDING 
+
+  fastjet::ClusterSequence clustSeq(inputParticles, plugin);
+
+//***************************** JETS EXTRACTION AND CORRECTION
+
+  // Here we extract inclusive jets with pt > ptmin, sorted by pt 
+  Double_t ptMin = header->GetMinJetPt(); 
+  vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
+  vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
+
+  for (Int_t j = 0; j < jets.size(); j++)
+  {
+    cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
+    cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<endl;
+    cout<<"px = "<<jets[j].px()<<endl;
+    cout<<"py = "<<jets[j].py()<<endl;
+    cout<<"pz = "<<jets[j].pz()<<endl;
+    cout<<"e = "<<jets[j].E()<<endl;
+    cout<<"******************"<<endl;
+    cout<<"******************"<<endl;
+    cout<<"******************"<<endl;
+
+    // Go to write AOD info
+    AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
+    if(debug) aodjet.Print("");
+    AddJet(aodjet);
+  }
+}
+//____________________________________________________________________________
+
+void AliSISConeJetFinder::WriteJHeaderToFile()
+{
+  fHeader->Write();
+}
diff --git a/JETAN/AliSISConeJetFinder.h b/JETAN/AliSISConeJetFinder.h
new file mode 100644 (file)
index 0000000..974b6d7
--- /dev/null
@@ -0,0 +1,58 @@
+#ifndef ALISISCONEJETFINDER_H
+#define ALISISCONEJETFINDER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+
+//---------------------------------------------------------------------
+// SISCone (FastJet v2.3.4) finder algorithm interface
+//
+// Author: swensy.jangal@ires.in2p3.fr
+//  
+//---------------------------------------------------------------------
+
+// FastJet classes 
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/PseudoJet.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 "AliFastJetHeader.h"
+#include "AliJetFinder.h"
+
+
+using namespace std;
+
+class AliSISConeJetFinder : public AliJetFinder
+{
+ public:
+
+  AliSISConeJetFinder();
+  ~AliSISConeJetFinder();
+
+  void FindJets(); 
+
+  // others
+  void WriteJHeaderToFile();
+  
+  protected:
+  AliSISConeJetFinder(const AliSISConeJetFinder& rfj);
+  AliSISConeJetFinder& operator = (const AliSISConeJetFinder& rsfj);
+
+  ClassDef(AliSISConeJetFinder,2)
+};
+
+#endif
diff --git a/JETAN/AliSISConeJetHeader.cxx b/JETAN/AliSISConeJetHeader.cxx
new file mode 100644 (file)
index 0000000..6e3ab66
--- /dev/null
@@ -0,0 +1,84 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+
+//---------------------------------------------------------------------
+// SISCone (FastJet v2.3.4) finder algorithm interface
+// Finder Header Class 
+// Author: swensy.jangal@ires.in2p3.fr
+//---------------------------------------------------------------------
+
+#include <Riostream.h>
+#include <TMath.h>
+
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/JetDefinition.hh"
+
+#include "AliSISConeJetHeader.h"
+
+ClassImp(AliSISConeJetHeader)
+
+////////////////////////////////////////////////////////////////////////
+
+AliSISConeJetHeader::AliSISConeJetHeader():
+    AliJetHeader("AliSISConeJetHeader"),
+    fActiveAreaRepeats(1),
+    fCaching(0),
+    fConeRadius(0.4),
+    fEffectiveRFact(1),
+    fGhostArea(0.05),
+    fGhostEtaMax(2.0),
+    fMinJetPt(0),
+    fNPassMax(0),
+    fOverlapThreshold(0.5),
+    fPhiMax(TMath::TwoPi()),
+    fPhiMin(0),
+    fRapMax(-0.9),
+    fRapMin(0.9),
+    fPtProtoJetMin(2),
+    fSplitMergeScaleNumber(0),
+    fSplitMergeStoppingScale(0)
+
+{
+  // Constructor
+}
+
+////////////////////////////////////////////////////////////////////////
+
+void AliSISConeJetHeader::PrintParameters() const
+{
+  // prints out parameters of jet algorithm
+
+  cout << "SISConeJet algorithm  parameters:"<<endl;
+
+  cout<<"Cone Radius = "<<fConeRadius<<endl;
+  cout<<"Overlap parameter = "<<fOverlapThreshold<<endl;
+  cout<<"Maximum number of runs = "<<fNPassMax<<endl;
+  cout<<"Pt min of protojets  = "<<fPtProtoJetMin<<endl;
+  cout<<"Do we record cones of these events ? (0 = no, 1 = yes) = "<<fCaching<<endl;
+
+  cout << "Background subtraction parameters :" <<endl;
+  //cout<<"Kind of area used = "<<<<endl;
+  cout<<"Eta max in which ghosts wil be generated = "<<fGhostEtaMax<<endl;
+  cout<<"Ghost area = "<<fGhostArea<<endl;
+  cout<<"Background will be studied in ["<<fRapMin<<","<<fRapMax<<"] in eta and ["<<fPhiMin<<","<<fPhiMax<<"] in phi"<<endl;
+  //cout<<"Kind of recombination for split/merge procedure = "<<<<endl;
+  //cout<<"Stopping scale for split/merge procedure = "<<<<endl;
+  cout<<"Do we repeat active area calculus? (0 = no, 1 = yes) = "<<fActiveAreaRepeats<<endl;
+
+  cout<<"Jets PtMin  = "<<fMinJetPt<<endl;
+  
+}
diff --git a/JETAN/AliSISConeJetHeader.h b/JETAN/AliSISConeJetHeader.h
new file mode 100644 (file)
index 0000000..ff24a76
--- /dev/null
@@ -0,0 +1,91 @@
+#ifndef ALISISCONEJETHEADER_H
+#define ALISISCONEJETHEADER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+//---------------------------------------------------------------------
+// SISCone (FastJet v2.3.4) finder algorithm interface
+// Finder Header Class 
+// Author: swensy.jangal@ires.in2p3.fr
+//---------------------------------------------------------------------
+#include "fastjet/AreaDefinition.hh"
+#include "fastjet/ClusterSequenceArea.hh"
+#include "fastjet/JetDefinition.hh"
+#include "AliJetHeader.h"
+
+class AliSISConeJetHeader : public AliJetHeader
+{
+ public:
+  AliSISConeJetHeader();
+  virtual ~AliSISConeJetHeader() { }
+
+  // Getters
+  Int_t                        GetActiveAreaRepeats()          {return fActiveAreaRepeats;}
+  Int_t                        GetAreaTypeNumber()             {return fAreaTypeNumber;}
+  Int_t                        GetNPassMax()                   {return fNPassMax;}
+  Int_t                        GetSplitMergeScale()            {return fSplitMergeScaleNumber;}
+
+  Double_t                     GetGhostEtaMax()                {return fGhostEtaMax;}
+  Double_t                     GetGhostArea()                  {return fGhostArea;}
+  Double_t                     GetEffectiveRFact()             {return fEffectiveRFact;}
+  Double_t                     GetRapMax()                     {return fRapMax;}
+  Double_t                     GetRapMin()                     {return fRapMin;}
+  Double_t                     GetPhiMax()                     {return fPhiMax;}
+  Double_t                     GetPhiMin()                     {return fPhiMin;}
+  Double_t                     GetConeRadius()                 {return fConeRadius;}
+  Double_t                     GetOverlapThreshold()           {return fOverlapThreshold;}
+  Double_t                     GetPtProtojetMin()              {return fPtProtoJetMin;}
+  Double_t                     GetCaching()                    {return fCaching;}
+  Double_t                     GetSplitMergeStoppingScale()    {return fSplitMergeStoppingScale;}
+  Double_t                     GetMinJetPt()                   {return fMinJetPt;}  
+
+  // Setters
+  void SetCaching(Bool_t value)                        {fCaching = value;}
+  void SetComment(TString com)                         {fComment=com;}
+  void SetGhostEtaMax(Double_t f)                      {fGhostEtaMax = f;}
+  void SetGhostArea(Double_t f)                        {fGhostArea = f;}
+  void SetActiveAreaRepeats(Int_t f)                   {fActiveAreaRepeats =f;}
+  void SetAreaTypeNumber(Int_t f)                      {fAreaTypeNumber = f;}
+  void SetEffectiveRFact(Double_t value)               {fEffectiveRFact = value;}       
+  void SetConeRadius(Double_t value)                   {fConeRadius = value;}
+  void SetMinJetPt(Double_t value)                     {fMinJetPt = value;}
+  void SetNPassMax(Int_t value)                        {fNPassMax = value;}
+  void SetOverlapThreshold(Double_t value)             {fOverlapThreshold = value;}
+  void SetPhiRange(Double_t fmin, Double_t fmax)       {fPhiMin = fmin; fPhiMax = fmax;}
+  void SetptProtojetMin(Double_t value)                {fPtProtoJetMin = value;}
+  void SetRapRange(Double_t fmin, Double_t fmax)       {fRapMin = fmin; fRapMax = fmax;}
+  void SetSplitMergeScale(Int_t value)                 {fSplitMergeScaleNumber = value;}
+  void SetSplitMergeStoppingScale(Double_t value)      {fSplitMergeStoppingScale = value;}       
+
+
+  // others
+  void PrintParameters() const;
+
+ protected:
+
+  Bool_t fCaching;                    // Do we record found cones for this set of data?
+
+  Int_t fActiveAreaRepeats;           // How many times do you want to caculate active areas?
+  Int_t fAreaTypeNumber;              // Kind of area
+  Int_t fNPassMax;                    // maximum number of passes
+  Int_t fSplitMergeScaleNumber;       // Kind of recombination in split/merge procedure, there's only one
+  
+  Double_t fConeRadius;               // Cone radius
+  Double_t fEffectiveRFact;           // Radius for Voronoi diagram
+  Double_t fGhostEtaMax;              // Maximum eta in which a ghost can be generated
+  Double_t fGhostArea;                // Area of one ghost
+  Double_t fMinJetPt;                 // Jet minimum energy
+  Double_t fOverlapThreshold;         // overlap parameter
+  Double_t fPhiMax, fPhiMin;          // Phi range
+  Double_t fPtProtoJetMin;            // pT min of protojets
+  Double_t fRapMax, fRapMin;          // Eta range
+  Double_t fSplitMergeStoppingScale;  // Stopping scale for split/merge procedure in case of area calculus
+
+  ClassDef(AliSISConeJetHeader,2)
+};
+#endif
index 6654254..e380b4c 100644 (file)
@@ -49,6 +49,8 @@
 #ifdef WITHFASTJET
 #pragma        link C++ class AliFastJetFinder+;
 #pragma        link C++ class AliFastJetHeader+;
+#pragma        link C++ class AliSISConeJetFinder+;
+#pragma        link C++ class AliSISConeJetHeader+;
 #endif
 
 #endif
index ad38bb2..5e4b37b 100644 (file)
@@ -29,7 +29,8 @@ SRCS =        AliJet.cxx AliJetHeader.cxx \
 
 
 ifneq ($(FASTJET),)
-       SRCS+= AliFastJetFinder.cxx AliFastJetHeader.cxx 
+       SRCS+= AliFastJetFinder.cxx AliFastJetHeader.cxx \
+       AliSISConeJetFinder.cxx AliSISConeJetHeader.cxx
 endif