Update responsibles for MCH, MTR, HMP
[u/mrichter/AliRoot.git] / JETAN / AliSISConeJetFinder.cxx
CommitLineData
2551af9d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
139cbd96 16/* $Id$ */
2551af9d 17
18//---------------------------------------------------------------------
8838ab7a 19// FastJet v2.3.4 finder algorithm interface
2551af9d 20//
21// Author: swensy.jangal@ires.in2p3.fr
8838ab7a 22//
139cbd96 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)
2551af9d 25//---------------------------------------------------------------------
26
27#include <Riostream.h>
139cbd96 28
29#include "AliFastJetInput.h"
30#include "AliFastJetBkg.h"
31#include "AliAODJetEventBackground.h"
32#include "AliAODJet.h"
2551af9d 33#include "AliSISConeJetFinder.h"
34#include "AliSISConeJetHeader.h"
35
36#include "fastjet/AreaDefinition.hh"
37#include "fastjet/ClusterSequenceArea.hh"
38#include "fastjet/JetDefinition.hh"
39#include "fastjet/PseudoJet.hh"
8838ab7a 40// get info on how fastjet was configured
2551af9d 41#include "fastjet/config.h"
42
d2eb0695 43#if defined(ENABLE_PLUGIN_SISCONE) || defined(FASTJET_ENABLE_PLUGIN_SISCONE)
2551af9d 44#include "fastjet/SISConePlugin.hh"
45#endif
46
2551af9d 47#include<vector>
2551af9d 48
2551af9d 49
139cbd96 50using namespace std;
8838ab7a 51
2551af9d 52ClassImp(AliSISConeJetFinder)
53
139cbd96 54////////////////////////////////////////////////////////////////////////
2551af9d 55
56AliSISConeJetFinder::AliSISConeJetFinder():
139cbd96 57 AliJetFinder(),
58 fInputFJ(new AliFastJetInput()),
59 fJetBkg(new AliFastJetBkg())
2551af9d 60{
61 // Constructor
62}
63
64//____________________________________________________________________________
2551af9d 65AliSISConeJetFinder::~AliSISConeJetFinder()
66{
67 // destructor
139cbd96 68 delete fInputFJ;
69 delete fJetBkg;
70
2551af9d 71}
72
73//______________________________________________________________________________
74void AliSISConeJetFinder::FindJets()
75{
139cbd96 76 // run the SISCone Jet finder
77
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
82
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__);
87 return;
88 }
2551af9d 89
139cbd96 90 //------------------- SISCONE PLUGIN CONFIGURATION ----------------------
91 // Look for SISCone parameters in the header and definition of the plugin.
2551af9d 92
139cbd96 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?
98 // For bckg
99 double rBkgParam = header->GetRparamBkg();
100 fastjet::Strategy strategy = header->GetStrategy();
101 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
2551af9d 102
139cbd96 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
2551af9d 105
106 fastjet::JetDefinition::Plugin * plugin;
107 plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
108
139cbd96 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
8838ab7a 112
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.
119
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;
127
128 fastjet::AreaDefinition areaDef;
129
130 if (areaTypeNumber < 5)
73faae2f 131 {
139cbd96 132 fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
133 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
134 }
73faae2f 135
139cbd96 136 if (areaTypeNumber == 5)
73faae2f 137 {
139cbd96 138 Double_t effectiveRFact = header->GetEffectiveRFact();
139 fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
140 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
73faae2f 141 }
73faae2f 142
139cbd96 143 //------------------- JETS FINDING AND EXTRACTION ----------------------
144 fastjet::ClusterSequenceArea clust_seq(inputParticles, plugin, areaDef);
73faae2f 145
139cbd96 146 vector<fastjet::PseudoJet> jets;
8838ab7a 147
139cbd96 148 if (bgMode == 1)// BG subtraction mode
8838ab7a 149 {
139cbd96 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);
158
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);
170 if(debug>0){
171 cout << "--------------------------------------------------------" << endl;
172 cout << comment << endl;
173 cout << "--------------------------------------------------------" << endl;
174 }
8838ab7a 175
139cbd96 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);
8838ab7a 180
139cbd96 181 //------------------- BACKGROUND SUBTRACTION ----------------------
182
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);
2551af9d 189
139cbd96 190 // Extract rho and sigma
191 Double_t rho = 0.;
192 Double_t sigma = 0.;
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);
2551af9d 197
139cbd96 198 // subtract background and extract jets bkg subtracted
199 vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptMin);
f92ba4be 200
139cbd96 201 // sort jets into increasing pt
202 jets = sorted_by_pt(subJets);
8838ab7a 203
139cbd96 204 }
205 else // No BG subtraction
8838ab7a 206 {
139cbd96 207
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);
215 if(debug>0){
216 cout << "--------------------------------------------------------" << endl;
217 cout << comment << endl;
218 cout << "--------------------------------------------------------" << endl;
9b72221b 219 }
139cbd96 220 //header->PrintParameters();
221
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);
f92ba4be 226
8838ab7a 227 }
95345ec0 228
139cbd96 229 //------------------- JET AND TRACK STORAGE ----------------------
230 for (size_t j = 0; j < jets.size(); j++) { // loop for jets
231
232 double area = clust_seq.area(jets[j]);
233 double areaError = clust_seq.area_error(jets[j]);
234
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);
236
237 vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
238 int nCon= constituents.size();
239 TArrayI ind(nCon);
240
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
246
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("");
252
253 for (int i=0; i < nCon; i++)
254 {
255 fastjet::PseudoJet mPart=constituents[i];
256 ind[i]=mPart.user_index();
257
258 // Jet constituents (charged tracks) added to the AliAODJet
259 AliJetCalTrkEvent* calEvt = GetCalTrkEvent();
260 for(Int_t itrack=0; itrack<calEvt->GetNCalTrkTracks(); itrack++)
261 {
262 if(itrack==ind[i])
263 {
264 TObject *track = calEvt->GetCalTrkTrack(itrack)->GetTrackObject();
265 aodjet.AddTrack(track);
266 }
267 }
268 } // End loop on Constituents
269
270 AddJet(aodjet);
271
272 } // end loop for jets
273
95345ec0 274 delete plugin;
275
2551af9d 276}
277
278//____________________________________________________________________________
139cbd96 279void AliSISConeJetFinder::WriteJHeaderToFile() const
8838ab7a 280{
139cbd96 281 // Write Jet Header To File(
282 fHeader->Write();
8838ab7a 283}
284
285//____________________________________________________________________________
139cbd96 286Bool_t AliSISConeJetFinder::ProcessEvent()
8838ab7a 287{
139cbd96 288 // Process one event
289 // Charged only or charged+neutral jets
8838ab7a 290
139cbd96 291 fInputFJ->SetHeader(fHeader);
292 fInputFJ->SetCalTrkEvent(GetCalTrkEvent());
293 fInputFJ->FillInput();
294
295 // Jets
296 FindJets();
297
298 // Background
299 if( fAODEvBkg){
300 fJetBkg->SetHeader(fHeader);
301 Double_t sigma1 = 0,meanarea1= 0,sigma2 = 0,meanarea2 = 0;
302 Double_t bkg1 = 0,bkg2 = 0;
303
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);
309 }
8838ab7a 310
139cbd96 311 Reset();
312 return kTRUE;
8838ab7a 313
2551af9d 314}