]>
Commit | Line | Data |
---|---|---|
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 | 50 | using namespace std; |
8838ab7a | 51 | |
2551af9d | 52 | ClassImp(AliSISConeJetFinder) |
53 | ||
139cbd96 | 54 | //////////////////////////////////////////////////////////////////////// |
2551af9d | 55 | |
56 | AliSISConeJetFinder::AliSISConeJetFinder(): | |
139cbd96 | 57 | AliJetFinder(), |
58 | fInputFJ(new AliFastJetInput()), | |
59 | fJetBkg(new AliFastJetBkg()) | |
2551af9d | 60 | { |
61 | // Constructor | |
62 | } | |
63 | ||
64 | //____________________________________________________________________________ | |
2551af9d | 65 | AliSISConeJetFinder::~AliSISConeJetFinder() |
66 | { | |
67 | // destructor | |
139cbd96 | 68 | delete fInputFJ; |
69 | delete fJetBkg; | |
70 | ||
2551af9d | 71 | } |
72 | ||
73 | //______________________________________________________________________________ | |
74 | void 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 | ||
8838ab7a | 278 | //____________________________________________________________________________ |
139cbd96 | 279 | void AliSISConeJetFinder::WriteJHeaderToFile() const |
8838ab7a | 280 | { |
139cbd96 | 281 | // Write Jet Header To File( |
282 | fHeader->Write(); | |
8838ab7a | 283 | } |
284 | ||
285 | //____________________________________________________________________________ | |
139cbd96 | 286 | Bool_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 | } |