1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 //---------------------------------------------------------------------
19 // FastJet v2.3.4 finder algorithm interface
21 // Author: swensy.jangal@ires.in2p3.fr
23 // ** 2011 magali.estienne@subatech.in2p3.fr & alexandre.shabetai@cern.ch
24 // Modified accordingly to reader/finder splitting and new handling of neutral information (via FastJetInput)
25 //---------------------------------------------------------------------
27 #include <Riostream.h>
29 #include "AliFastJetInput.h"
30 #include "AliFastJetBkg.h"
31 #include "AliAODJetEventBackground.h"
32 #include "AliAODJet.h"
33 #include "AliSISConeJetFinder.h"
34 #include "AliSISConeJetHeader.h"
36 #include "fastjet/AreaDefinition.hh"
37 #include "fastjet/ClusterSequenceArea.hh"
38 #include "fastjet/JetDefinition.hh"
39 #include "fastjet/PseudoJet.hh"
40 // get info on how fastjet was configured
41 #include "fastjet/config.h"
43 #if defined(ENABLE_PLUGIN_SISCONE) || defined(FASTJET_ENABLE_PLUGIN_SISCONE)
44 #include "fastjet/SISConePlugin.hh"
52 ClassImp(AliSISConeJetFinder)
54 ////////////////////////////////////////////////////////////////////////
56 AliSISConeJetFinder::AliSISConeJetFinder():
58 fInputFJ(new AliFastJetInput()),
59 fJetBkg(new AliFastJetBkg())
64 //____________________________________________________________________________
65 AliSISConeJetFinder::~AliSISConeJetFinder()
73 //______________________________________________________________________________
74 void AliSISConeJetFinder::FindJets()
76 // run the SISCone Jet finder
78 // Pick up siscone header
79 AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader;
80 Int_t debug = header->GetDebug(); // debug option
81 Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
83 // Read input particles
84 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
85 if(inputParticles.size()==0){
86 if(debug>0) Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
90 //------------------- SISCONE PLUGIN CONFIGURATION ----------------------
91 // Look for SISCone parameters in the header and definition of the plugin.
93 Double_t coneRadius = header->GetConeRadius(); // cone radius
94 Double_t overlapThreshold = header->GetOverlapThreshold(); // overlap parameter
95 Int_t nPassMax = header->GetNPassMax(); // maximum number of passes
96 Double_t ptProtoJetMin = header->GetPtProtojetMin(); // pT min of protojets
97 Double_t caching = header->GetCaching(); // do we record found cones for this set of data?
99 double rBkgParam = header->GetRparamBkg();
100 fastjet::Strategy strategy = header->GetStrategy();
101 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
103 // if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
104 // Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets
106 fastjet::JetDefinition::Plugin * plugin;
107 plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
109 //------------------- CHOICE OF JET AREA ----------------------
110 // Definition of jet areas for background subtraction
111 // For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
113 Double_t ghostEtamax = header->GetGhostEtaMax(); // maximum eta in which a ghost can be generated
114 Double_t ghostArea = header->GetGhostArea(); // area of a ghost
115 Int_t activeAreaRepeats = header->GetActiveAreaRepeats(); // do we repeat area calculation?
116 Double_t gridScatter = header->GetGridScatter(); // fractional random fluctuations of the position of the ghosts on the y-phi grid
117 Double_t ktScatter = header->GetKtScatter(); // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid
118 Double_t meanGhostKt = header->GetMeanGhostKt(); // average transverse momentum of the ghosts.
120 Double_t areaTypeNumber = header->GetAreaTypeNumber(); // the number determines jet area type
121 fastjet::AreaType areaType = fastjet::active_area;
122 if (areaTypeNumber == 1) areaType = fastjet::active_area;
123 if (areaTypeNumber == 2) areaType = fastjet::active_area_explicit_ghosts;
124 if (areaTypeNumber == 3) areaType = fastjet::one_ghost_passive_area;
125 if (areaTypeNumber == 4) areaType = fastjet::passive_area;
126 if (areaTypeNumber == 5) areaType = fastjet::voronoi_area;
128 fastjet::AreaDefinition areaDef;
130 if (areaTypeNumber < 5)
132 fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
133 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
136 if (areaTypeNumber == 5)
138 Double_t effectiveRFact = header->GetEffectiveRFact();
139 fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
140 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
143 //------------------- JETS FINDING AND EXTRACTION ----------------------
144 fastjet::ClusterSequenceArea clust_seq(inputParticles, plugin, areaDef);
146 vector<fastjet::PseudoJet> jets;
148 if (bgMode == 1)// BG subtraction mode
150 //------------------- CLUSTER JETS FINDING FOR RHO ESTIMATION ----------------------
151 // run the jet clustering with the above jet definition
152 fastjet::JetAlgorithm algorithmBkg = fastjet::kt_algorithm;
153 Int_t algo = header->GetBGAlgorithm();
154 if (algo == 0) algorithmBkg = fastjet::kt_algorithm;
155 if (algo == 1) algorithmBkg = fastjet::cambridge_algorithm;
156 fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy);
157 fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef);
159 // save a comment in the header
160 TString comment = "Running Siscone algorithm with the following setup. ";
161 comment+= "Jet definition: ";
162 // comment+= TString(plugin.description());
163 comment+= "Jet bckg definition: ";
164 comment+= TString(jetDefBkg.description());
165 comment+= ". Area definition: ";
166 comment+= TString(areaDef.description());
167 comment+= ". Strategy adopted by FastJet and bkg: ";
168 comment+= TString(clust_seq.strategy_string());
169 header->SetComment(comment);
171 cout << "--------------------------------------------------------" << endl;
172 cout << comment << endl;
173 cout << "--------------------------------------------------------" << endl;
176 // Here we extract inclusive jets with pt > ptmin, sorted by pt
177 Double_t ptMin = header->GetMinJetPt();
178 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(); // ptMin removed
179 // vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
181 //------------------- BACKGROUND SUBTRACTION ----------------------
183 // Set the rapidity-azimuth range within which to study background
184 Double_t rapMin = header->GetRapMin();
185 Double_t rapMax = header->GetRapMax();
186 Double_t phiMin = header->GetPhiMin();
187 Double_t phiMax = header->GetPhiMax();
188 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
190 // Extract rho and sigma
193 Double_t meanarea = 0.;
194 Bool_t kUse4VectorArea = header->Use4VectorArea();
195 vector<fastjet::PseudoJet> bkgJets = clust_seq_bkg.inclusive_jets();
196 clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, kUse4VectorArea, rho, sigma, meanarea, false);
198 // subtract background and extract jets bkg subtracted
199 vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptMin);
201 // sort jets into increasing pt
202 jets = sorted_by_pt(subJets);
205 else // No BG subtraction
208 // save a comment in the header
209 TString comment = "Running Siscone algorithm with the following setup. ";
210 comment+= "Jet definition: ";
211 // comment+= TString(plugin.description());
212 comment+= ". Strategy adopted by FastJet: ";
213 comment+= TString(clust_seq.strategy_string());
214 header->SetComment(comment);
216 cout << "--------------------------------------------------------" << endl;
217 cout << comment << endl;
218 cout << "--------------------------------------------------------" << endl;
220 //header->PrintParameters();
222 // Here we extract inclusive jets with pt > ptmin, sorted by pt
223 Double_t ptMin = header->GetMinJetPt();
224 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptMin);
225 jets = sorted_by_pt(inclusiveJets);
229 //------------------- JET AND TRACK STORAGE ----------------------
230 for (size_t j = 0; j < jets.size(); j++) { // loop for jets
232 double area = clust_seq.area(jets[j]);
233 double areaError = clust_seq.area_error(jets[j]);
235 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);
237 vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
238 int nCon= constituents.size();
241 if ((jets[j].eta() > (header->GetJetEtaMax())) ||
242 (jets[j].eta() < (header->GetJetEtaMin())) ||
243 (jets[j].phi() > (header->GetJetPhiMax())) ||
244 (jets[j].phi() < (header->GetJetPhiMin())) ||
245 (jets[j].perp() < header->GetMinJetPt())) continue; // acceptance eta range and etmin
247 // go to write AOD info
248 AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
249 aodjet.SetEffArea(area,areaError);
250 //cout << "Printing jet " << endl;
251 if(debug>0) aodjet.Print("");
253 for (int i=0; i < nCon; i++)
255 fastjet::PseudoJet mPart=constituents[i];
256 ind[i]=mPart.user_index();
258 // Jet constituents (charged tracks) added to the AliAODJet
259 AliJetCalTrkEvent* calEvt = GetCalTrkEvent();
260 for(Int_t itrack=0; itrack<calEvt->GetNCalTrkTracks(); itrack++)
264 TObject *track = calEvt->GetCalTrkTrack(itrack)->GetTrackObject();
265 aodjet.AddTrack(track);
268 } // End loop on Constituents
272 } // end loop for jets
278 //____________________________________________________________________________
279 void AliSISConeJetFinder::WriteJHeaderToFile() const
281 // Write Jet Header To File(
285 //____________________________________________________________________________
286 Bool_t AliSISConeJetFinder::ProcessEvent()
289 // Charged only or charged+neutral jets
291 fInputFJ->SetHeader(fHeader);
292 fInputFJ->SetCalTrkEvent(GetCalTrkEvent());
293 fInputFJ->FillInput();
300 fJetBkg->SetHeader(fHeader);
301 Double_t sigma1 = 0,meanarea1= 0,sigma2 = 0,meanarea2 = 0;
302 Double_t bkg1 = 0,bkg2 = 0;
304 fJetBkg->SetFastJetInput(fInputFJ);
305 fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
306 fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
307 fAODEvBkg->SetBackground(0,bkg1,sigma1,meanarea1);
308 fAODEvBkg->SetBackground(1,bkg2,sigma2,meanarea2);