]>
Commit | Line | Data |
---|---|---|
d89b8229 | 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 | ||
16 | /* $Id$ */ | |
17 | ||
18 | //--------------------------------------------------------------------- | |
19 | // FastJet v2.3.4 finder algorithm interface | |
20 | // | |
21 | // Author: swensy.jangal@ires.in2p3.fr | |
22 | // | |
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 | //--------------------------------------------------------------------- | |
26 | ||
27 | #include <Riostream.h> | |
28 | ||
29 | #include "AliFastJetInput.h" | |
30 | #include "AliFastJetBkg.h" | |
31 | #include "AliAODJetEventBackground.h" | |
32 | #include "AliAODJet.h" | |
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" | |
40 | // get info on how fastjet was configured | |
41 | #include "fastjet/config.h" | |
42 | ||
43 | #ifdef ENABLE_PLUGIN_SISCONE | |
44 | #include "fastjet/SISConePlugin.hh" | |
45 | #endif | |
46 | ||
47 | #include<vector> | |
48 | ||
49 | ||
50 | using namespace std; | |
51 | ||
52 | ClassImp(AliSISConeJetFinder) | |
53 | ||
54 | //////////////////////////////////////////////////////////////////////// | |
55 | ||
56 | AliSISConeJetFinder::AliSISConeJetFinder(): | |
57 | AliJetFinder(), | |
58 | fInputFJ(new AliFastJetInput()), | |
59 | fJetBkg(new AliFastJetBkg()) | |
60 | { | |
61 | // Constructor | |
62 | } | |
63 | ||
64 | //____________________________________________________________________________ | |
65 | AliSISConeJetFinder::~AliSISConeJetFinder() | |
66 | { | |
67 | // destructor | |
68 | delete fInputFJ; | |
69 | delete fJetBkg; | |
70 | ||
71 | } | |
72 | ||
73 | //______________________________________________________________________________ | |
74 | void AliSISConeJetFinder::FindJets() | |
75 | { | |
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 | } | |
89 | ||
90 | //------------------- SISCONE PLUGIN CONFIGURATION ---------------------- | |
91 | // Look for SISCone parameters in the header and definition of the plugin. | |
92 | ||
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(); | |
102 | ||
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 | |
105 | ||
106 | fastjet::JetDefinition::Plugin * plugin; | |
107 | plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching); | |
108 | ||
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 | |
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) | |
131 | { | |
132 | fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt); | |
133 | areaDef = fastjet::AreaDefinition(areaType,ghostSpec); | |
134 | } | |
135 | ||
136 | if (areaTypeNumber == 5) | |
137 | { | |
138 | Double_t effectiveRFact = header->GetEffectiveRFact(); | |
139 | fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact); | |
140 | areaDef = fastjet::AreaDefinition(areaType,ghostSpec); | |
141 | } | |
142 | ||
143 | //------------------- JETS FINDING AND EXTRACTION ---------------------- | |
144 | fastjet::ClusterSequenceArea clust_seq(inputParticles, plugin, areaDef); | |
145 | ||
146 | vector<fastjet::PseudoJet> jets; | |
147 | ||
148 | if (bgMode == 1)// BG subtraction mode | |
149 | { | |
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 | } | |
175 | ||
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); | |
180 | ||
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); | |
189 | ||
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); | |
197 | ||
198 | // subtract background and extract jets bkg subtracted | |
199 | vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptMin); | |
200 | ||
201 | // sort jets into increasing pt | |
202 | jets = sorted_by_pt(subJets); | |
203 | ||
204 | } | |
205 | else // No BG subtraction | |
206 | { | |
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; | |
219 | } | |
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); | |
226 | ||
227 | } | |
228 | ||
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 | ||
274 | delete plugin; | |
275 | ||
276 | } | |
277 | ||
278 | //____________________________________________________________________________ | |
279 | void AliSISConeJetFinder::WriteJHeaderToFile() const | |
280 | { | |
281 | // Write Jet Header To File( | |
282 | fHeader->Write(); | |
283 | } | |
284 | ||
285 | //____________________________________________________________________________ | |
286 | Bool_t AliSISConeJetFinder::ProcessEvent() | |
287 | { | |
288 | // Process one event | |
289 | // Charged only or charged+neutral jets | |
290 | ||
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 | } | |
310 | ||
311 | Reset(); | |
312 | return kTRUE; | |
313 | ||
314 | } |