* provided "as is" without express or implied warranty. *
**************************************************************************/
+/* $Id$ */
//---------------------------------------------------------------------
// FastJet v2.3.4 finder algorithm interface
//
// Author: swensy.jangal@ires.in2p3.fr
//
-// Last modification: Neutral cell energy included in the jet reconstruction
-//
-// Author: Magali.estienne@subatech.in2p3.fr
+// ** 2011 magali.estienne@subatech.in2p3.fr & alexandre.shabetai@cern.ch
+// Modified accordingly to reader/finder splitting and new handling of neutral information (via FastJetInput)
//---------------------------------------------------------------------
-
#include <Riostream.h>
-#include <TArrayF.h>
-#include <TFile.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TLorentzVector.h>
-#include <TRandom.h>
-#include <TClonesArray.h>
-
-#include "AliHeader.h"
-#include "AliJetKineReader.h"
-#include "AliJetReader.h"
-#include "AliJetReaderHeader.h"
-#include "AliJetUnitArray.h"
+
+#include "AliFastJetInput.h"
+#include "AliFastJetBkg.h"
+#include "AliAODJetEventBackground.h"
+#include "AliAODJet.h"
#include "AliSISConeJetFinder.h"
#include "AliSISConeJetHeader.h"
#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
+#if defined(ENABLE_PLUGIN_SISCONE) || defined(FASTJET_ENABLE_PLUGIN_SISCONE)
#include "fastjet/SISConePlugin.hh"
#endif
-#include<sstream> // needed for internal io
#include<vector>
-#include <cmath>
-using namespace std;
+using namespace std;
ClassImp(AliSISConeJetFinder)
-
-//____________________________________________________________________________
+////////////////////////////////////////////////////////////////////////
AliSISConeJetFinder::AliSISConeJetFinder():
-AliJetFinder()
+ AliJetFinder(),
+ fInputFJ(new AliFastJetInput()),
+ fJetBkg(new AliFastJetBkg())
{
// Constructor
}
//____________________________________________________________________________
-
AliSISConeJetFinder::~AliSISConeJetFinder()
{
// destructor
+ delete fInputFJ;
+ delete fJetBkg;
+
}
//______________________________________________________________________________
void AliSISConeJetFinder::FindJets()
{
+ // run the SISCone Jet finder
- // Pick up siscone header
+ // Pick up siscone header
AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader;
- Int_t debug = header->GetDebug(); // debug option
- Int_t fOpt = fReader->GetReaderHeader()->GetDetector();
-
- // Check if we are reading AOD jets
- Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
+ Int_t debug = header->GetDebug(); // debug option
+ Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
+
+ // Read input particles
+ vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
+ if(inputParticles.size()==0){
+ if(debug>0) Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
+ return;
+ }
-//******************************** SISCONE PLUGIN CONFIGURATION
-// Here we look for SISCone parameters in the header and we define our plugin.
+ //------------------- SISCONE PLUGIN CONFIGURATION ----------------------
+ // Look for SISCone parameters in the header and definition of the 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?
+ 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?
+ // For bckg
+ double rBkgParam = header->GetRparamBkg();
+ fastjet::Strategy strategy = header->GetStrategy();
+ fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
-// if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
-// Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets
+ // if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
+ // Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets
fastjet::JetDefinition::Plugin * plugin;
plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
- vector<fastjet::PseudoJet> inputParticles;
-//******************************** 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.
-
- if(fOpt==0)
- {
- 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;
- delete plugin;
- return;
- }
-
- if(nIn == 0)// nIn = Number of particles in the event
- {
- if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
- delete plugin;
- return;
- }
-
- Float_t px,py,pz,en;
-
- // 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);
-
- }
- }
- else {
- TClonesArray* fUnit = fReader->GetUnitArray();
- if(fUnit == 0) { cout << "Could not get the momentum array" << endl; delete plugin; return; }
- Int_t nIn = fUnit->GetEntries();
- if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; delete plugin; return; }
- // Information extracted from fUnitArray
- // load input vectors and calculate total energy in array
- Float_t pt,eta,phi,theta,px,py,pz,en;
- Int_t ipart = 0;
- for(Int_t i=0; i<nIn; i++)
- {
- AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
-
- if(uArray->GetUnitEnergy()>0.){
-
- // It is not necessary anymore to cut on particle pt
- pt = uArray->GetUnitEnergy();
- eta = uArray->GetUnitEta();
- phi = uArray->GetUnitPhi();
- theta = EtaToTheta(eta);
- en = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta));
- px = TMath::Cos(phi)*pt;
- py = TMath::Sin(phi)*pt;
- pz = en*TMath::TanH(eta);
- if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
-
- fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object
- inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm
- inputParticles.push_back(inputPart); // back of the input_particles vector
- ipart++;
- }
- } // End loop on UnitArray
- }
-
-//******************************** CHOICE OF JET AREA
-// Here we determine jets area for subtracting background later
-// For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
+ //------------------- CHOICE OF JET AREA ----------------------
+ // Definition of jet areas for background subtraction
+ // For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
Double_t ghostEtamax = header->GetGhostEtaMax(); // maximum eta in which a ghost can be generated
Double_t ghostArea = header->GetGhostArea(); // area of a ghost
fastjet::AreaDefinition areaDef;
if (areaTypeNumber < 5)
- {
- fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
- areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
- }
+ {
+ fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
+ areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
+ }
if (areaTypeNumber == 5)
- {
- Double_t effectiveRFact = header->GetEffectiveRFact();
- fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
- areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
- }
-
-//********************************
-//********************************
-
- Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
+ {
+ Double_t effectiveRFact = header->GetEffectiveRFact();
+ fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
+ areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
+ }
-//********************************
-//********************************
+ //------------------- JETS FINDING AND EXTRACTION ----------------------
+ fastjet::ClusterSequenceArea clust_seq(inputParticles, plugin, areaDef);
-//********************************
-//******************************** BG SUBTRACTION
- if (bgMode == 1)// BG subtraction
- {
- //******************************** JETS FINDING AND EXTRACTION
- fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
- // 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);
-
- //***************************** BACKGROUND SUBTRACTION
-
- // Set the rapidity-azimuth range within which to study background
- Double_t rapMin = header->GetRapMin();
- Double_t rapMax = header->GetRapMax();
- Double_t phiMin = header->GetPhiMin();
- Double_t phiMax = header->GetPhiMax();
- fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
-
- // rho calculus from Cambridge/Aachen clustering (not from SISCONE as it gives too small areas)
- Int_t algo = header->GetBGAlgorithm();
- fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
- if (algo == 0) algorithm = fastjet::kt_algorithm;
- if (algo == 1) algorithm = fastjet::cambridge_algorithm;
- Double_t RRho = header->GetRForRho();
- fastjet::JetDefinition jetDefForRho(algorithm, RRho);
- fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefForRho, areaDef);
- Double_t rho = csForRho.median_pt_per_unit_area_4vector(range);
- cout<<"rho = "<<rho<<endl;
-
- // Vector of corrected jets
- vector<fastjet::PseudoJet> corrJets = clustSeq.subtracted_jets(rho, ptMin);
-
- //***************************** JETS DISPLAY
+ vector<fastjet::PseudoJet> jets;
- for (size_t j = 0; j < jets.size(); j++)
+ if (bgMode == 1)// BG subtraction mode
{
- // If the jet is only made of ghosts, continue.
- if (clustSeq.is_pure_ghost(jets[j]) == 1) continue;
+ //------------------- CLUSTER JETS FINDING FOR RHO ESTIMATION ----------------------
+ // run the jet clustering with the above jet definition
+ fastjet::JetAlgorithm algorithmBkg = fastjet::kt_algorithm;
+ Int_t algo = header->GetBGAlgorithm();
+ if (algo == 0) algorithmBkg = fastjet::kt_algorithm;
+ if (algo == 1) algorithmBkg = fastjet::cambridge_algorithm;
+ fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy);
+ fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef);
+
+ // save a comment in the header
+ TString comment = "Running Siscone algorithm with the following setup. ";
+ comment+= "Jet definition: ";
+ // comment+= TString(plugin.description());
+ comment+= "Jet bckg definition: ";
+ comment+= TString(jetDefBkg.description());
+ comment+= ". Area definition: ";
+ comment+= TString(areaDef.description());
+ comment+= ". Strategy adopted by FastJet and bkg: ";
+ comment+= TString(clust_seq.strategy_string());
+ header->SetComment(comment);
+ if(debug>0){
+ cout << "--------------------------------------------------------" << endl;
+ cout << comment << endl;
+ cout << "--------------------------------------------------------" << endl;
+ }
- // If the correction is > jet energy px = py = pz = e = 0
- if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue;
+ // Here we extract inclusive jets with pt > ptmin, sorted by pt
+ Double_t ptMin = header->GetMinJetPt();
+ vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(); // ptMin removed
+ // vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
- if(debug){
- cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
- cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<" area : "<<clustSeq.area(jets[j])<<endl;
- cout<<"e = "<<jets[j].E()<<endl;
-
- cout<<"********************************** Corrected jet(s)"<<endl;
- cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<corrJets[j].rap()<<" Phi : "<<corrJets[j].phi()<<" pT : "<<corrJets[j].perp()<<endl;
- cout<<"e = "<<corrJets[j].E()<<endl;
- }
- // Go to write AOD info
- AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E());
- if(debug) aodjet.Print("");
- AddJet(aodjet);
- }
- }
+ //------------------- BACKGROUND SUBTRACTION ----------------------
-//********************************
-//******************************** NO BG SUBTRACTION
+ // Set the rapidity-azimuth range within which to study background
+ Double_t rapMin = header->GetRapMin();
+ Double_t rapMax = header->GetRapMax();
+ Double_t phiMin = header->GetPhiMin();
+ Double_t phiMax = header->GetPhiMax();
+ fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
- if (bgMode == 0)// No BG subtraction
- {
- //******************************** JETS FINDING AND EXTRACTION
- fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
- // 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);
+ // Extract rho and sigma
+ Double_t rho = 0.;
+ Double_t sigma = 0.;
+ Double_t meanarea = 0.;
+ Bool_t kUse4VectorArea = header->Use4VectorArea();
+ vector<fastjet::PseudoJet> bkgJets = clust_seq_bkg.inclusive_jets();
+ clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, kUse4VectorArea, rho, sigma, meanarea, false);
- //***************************** JETS DISPLAY
+ // subtract background and extract jets bkg subtracted
+ vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptMin);
- for (size_t k = 0; k < jets.size(); k++)
+ // sort jets into increasing pt
+ jets = sorted_by_pt(subJets);
+
+ }
+ else // No BG subtraction
{
- if(debug)
- {
- cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
- cout<<"Jet number "<<k+1<<" : "<<"Rapidity : "<<jets[k].rap()<<" Phi : "<<jets[k].phi()<<" pT : "<<jets[k].perp()<<endl;
- cout<<"e = "<<jets[k].E()<<endl;
+
+ // save a comment in the header
+ TString comment = "Running Siscone algorithm with the following setup. ";
+ comment+= "Jet definition: ";
+ // comment+= TString(plugin.description());
+ comment+= ". Strategy adopted by FastJet: ";
+ comment+= TString(clust_seq.strategy_string());
+ header->SetComment(comment);
+ if(debug>0){
+ cout << "--------------------------------------------------------" << endl;
+ cout << comment << endl;
+ cout << "--------------------------------------------------------" << endl;
}
+ //header->PrintParameters();
+
+ // Here we extract inclusive jets with pt > ptmin, sorted by pt
+ Double_t ptMin = header->GetMinJetPt();
+ vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptMin);
+ jets = sorted_by_pt(inclusiveJets);
- // Go to write AOD info
- Double_t area = clustSeq.area(jets[k]);
- AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
- aodjet.SetEffArea(area,0);
- if(debug) aodjet.Print("");
- AddJet(aodjet);
}
- }
+ //------------------- JET AND TRACK STORAGE ----------------------
+ for (size_t j = 0; j < jets.size(); j++) { // loop for jets
+
+ double area = clust_seq.area(jets[j]);
+ double areaError = clust_seq.area_error(jets[j]);
+
+ if(debug>0) printf("Jet found %5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, areaError);
+
+ vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
+ int nCon= constituents.size();
+ TArrayI ind(nCon);
+
+ if ((jets[j].eta() > (header->GetJetEtaMax())) ||
+ (jets[j].eta() < (header->GetJetEtaMin())) ||
+ (jets[j].phi() > (header->GetJetPhiMax())) ||
+ (jets[j].phi() < (header->GetJetPhiMin())) ||
+ (jets[j].perp() < header->GetMinJetPt())) continue; // acceptance eta range and etmin
+
+ // go to write AOD info
+ AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
+ aodjet.SetEffArea(area,areaError);
+ //cout << "Printing jet " << endl;
+ if(debug>0) aodjet.Print("");
+
+ for (int i=0; i < nCon; i++)
+ {
+ fastjet::PseudoJet mPart=constituents[i];
+ ind[i]=mPart.user_index();
+
+ // Jet constituents (charged tracks) added to the AliAODJet
+ AliJetCalTrkEvent* calEvt = GetCalTrkEvent();
+ for(Int_t itrack=0; itrack<calEvt->GetNCalTrkTracks(); itrack++)
+ {
+ if(itrack==ind[i])
+ {
+ TObject *track = calEvt->GetCalTrkTrack(itrack)->GetTrackObject();
+ aodjet.AddTrack(track);
+ }
+ }
+ } // End loop on Constituents
+
+ AddJet(aodjet);
+
+ } // end loop for jets
+
delete plugin;
}
//____________________________________________________________________________
-
-Float_t AliSISConeJetFinder::EtaToTheta(Float_t arg)
+void AliSISConeJetFinder::WriteJHeaderToFile() const
{
- // return (180./TMath::Pi())*2.*atan(exp(-arg));
- return 2.*atan(exp(-arg));
-
-
+ // Write Jet Header To File(
+ fHeader->Write();
}
//____________________________________________________________________________
-
-void AliSISConeJetFinder::InitTask(TChain *tree)
+Bool_t AliSISConeJetFinder::ProcessEvent()
{
+ // Process one event
+ // Charged only or charged+neutral jets
- printf("SISCone jet finder initialization ******************");
- fReader->CreateTasks(tree);
-
-}
+ fInputFJ->SetHeader(fHeader);
+ fInputFJ->SetCalTrkEvent(GetCalTrkEvent());
+ fInputFJ->FillInput();
+
+ // Jets
+ FindJets();
+
+ // Background
+ if( fAODEvBkg){
+ fJetBkg->SetHeader(fHeader);
+ Double_t sigma1 = 0,meanarea1= 0,sigma2 = 0,meanarea2 = 0;
+ Double_t bkg1 = 0,bkg2 = 0;
+
+ fJetBkg->SetFastJetInput(fInputFJ);
+ fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
+ fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
+ fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+ fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2);
+ }
+ Reset();
+ return kTRUE;
-//____________________________________________________________________________
-
-void AliSISConeJetFinder::WriteJHeaderToFile() const
-{
- fHeader->Write();
}
-
-//____________________________________________________________________________
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff