]>
Commit | Line | Data |
---|---|---|
a17e6965 | 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 | ||
392c9b47 | 16 | |
a17e6965 | 17 | //--------------------------------------------------------------------- |
392c9b47 | 18 | // FastJet v2.3.4 finder algorithm interface |
8838ab7a | 19 | // Last modification: Neutral cell energy included in the jet reconstruction |
20 | // | |
21 | // Authors: Rafael.Diaz.Valdes@cern.ch | |
22 | // Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option) | |
392c9b47 | 23 | // |
a17e6965 | 24 | //--------------------------------------------------------------------- |
25 | ||
8838ab7a | 26 | |
a17e6965 | 27 | #include <Riostream.h> |
28 | #include <TLorentzVector.h> | |
29 | #include <TFile.h> | |
30 | #include <TH1F.h> | |
31 | #include <TH2F.h> | |
32 | #include <TArrayF.h> | |
33 | #include <TRandom.h> | |
34 | #include <TClonesArray.h> | |
35 | ||
36 | #include "AliFastJetFinder.h" | |
8838ab7a | 37 | #include "AliFastJetHeaderV1.h" |
a17e6965 | 38 | #include "AliJetReaderHeader.h" |
39 | #include "AliJetReader.h" | |
40 | #include "AliJet.h" | |
8838ab7a | 41 | #include "AliJetUnitArray.h" |
392c9b47 | 42 | |
43 | #include "fastjet/PseudoJet.hh" | |
44 | #include "fastjet/ClusterSequenceArea.hh" | |
45 | #include "fastjet/AreaDefinition.hh" | |
46 | #include "fastjet/JetDefinition.hh" | |
47 | // get info on how fastjet was configured | |
48 | #include "fastjet/config.h" | |
49 | ||
50 | #ifdef ENABLE_PLUGIN_SISCONE | |
51 | #include "fastjet/SISConePlugin.hh" | |
52 | #endif | |
53 | ||
54 | #include<sstream> // needed for internal io | |
55 | #include<vector> | |
56 | #include <cmath> | |
57 | ||
58 | using namespace std; | |
a17e6965 | 59 | |
60 | ||
9157b9c9 | 61 | ClassImp(AliFastJetFinder) |
8838ab7a | 62 | |
63 | ||
392c9b47 | 64 | //____________________________________________________________________________ |
a17e6965 | 65 | |
66 | AliFastJetFinder::AliFastJetFinder(): | |
392c9b47 | 67 | AliJetFinder() |
a17e6965 | 68 | { |
69 | // Constructor | |
70 | } | |
71 | ||
392c9b47 | 72 | //____________________________________________________________________________ |
a17e6965 | 73 | |
74 | AliFastJetFinder::~AliFastJetFinder() | |
a17e6965 | 75 | { |
76 | // destructor | |
77 | } | |
78 | ||
392c9b47 | 79 | //______________________________________________________________________________ |
a17e6965 | 80 | void AliFastJetFinder::FindJets() |
a17e6965 | 81 | { |
392c9b47 | 82 | |
392c9b47 | 83 | //pick up fastjet header |
8838ab7a | 84 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; |
85 | Bool_t debug = header->GetDebug(); // debug option | |
86 | Bool_t bgMode = header->GetBGMode(); // choose to subtract BG or not | |
87 | Int_t fOpt = fReader->GetReaderHeader()->GetDetector(); | |
392c9b47 | 88 | |
89 | // check if we are reading AOD jets | |
90 | TRefArray *refs = 0; | |
91 | Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader"); | |
92 | if (fromAod) { refs = fReader->GetReferences(); } | |
93 | ||
94 | // RUN ALGORITHM | |
95 | // read input particles ----------------------------- | |
96 | vector<fastjet::PseudoJet> input_particles; | |
8838ab7a | 97 | if(fOpt==0) |
98 | { | |
99 | TClonesArray *lvArray = fReader->GetMomentumArray(); | |
100 | if(lvArray == 0) { cout << "Could not get the momentum array" << endl; return; } | |
101 | Int_t nIn = lvArray->GetEntries(); | |
102 | if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; } | |
103 | fJets->SetNinput(nIn) ; // number of input objects | |
104 | Float_t px,py,pz,en; | |
105 | // load input vectors | |
106 | for(Int_t i = 0; i < nIn; i++){ // loop for all input particles | |
107 | TLorentzVector *lv = (TLorentzVector*) lvArray->At(i); | |
108 | px = lv->Px(); | |
109 | py = lv->Py(); | |
110 | pz = lv->Pz(); | |
111 | en = lv->Energy(); | |
112 | ||
113 | fastjet::PseudoJet input_part(px,py,pz,en); // create PseudoJet object | |
114 | input_part.set_user_index(i); //label the particle into Fastjet algortihm | |
115 | input_particles.push_back(input_part); // back of the input_particles vector | |
116 | } // end loop | |
117 | } | |
118 | else { | |
119 | TClonesArray* fUnit = fReader->GetUnitArray(); | |
120 | if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; } | |
121 | Int_t nCandidate = fReader->GetNumCandidate(); | |
122 | Int_t nIn = fUnit->GetEntries(); | |
123 | if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; } | |
124 | fJets->SetNinput(nCandidate); // number of input objects // ME | |
125 | // Information extracted from fUnitArray | |
126 | // load input vectors and calculate total energy in array | |
127 | Float_t pt,eta,phi,theta,px,py,pz,en; | |
128 | Int_t ipart = 0; | |
129 | for(Int_t i=0; i<nIn; i++) | |
130 | { | |
131 | AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i); | |
132 | ||
133 | if(uArray->GetUnitEnergy()>0.){ | |
134 | ||
135 | // It is not necessary anymore to cut on particle pt | |
136 | pt = uArray->GetUnitEnergy(); | |
137 | eta = uArray->GetUnitEta(); | |
138 | phi = uArray->GetUnitPhi(); | |
139 | theta = EtaToTheta(eta); | |
140 | en = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta)); | |
141 | px = TMath::Cos(phi)*pt; | |
142 | py = TMath::Sin(phi)*pt; | |
143 | pz = en*TMath::TanH(eta); | |
144 | if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl; | |
145 | ||
146 | fastjet::PseudoJet input_part(px,py,pz,en); // create PseudoJet object | |
147 | input_part.set_user_index(ipart); //label the particle into Fastjet algortihm | |
148 | input_particles.push_back(input_part); // back of the input_particles vector | |
149 | ipart++; | |
150 | } | |
151 | } // End loop on UnitArray | |
152 | } | |
392c9b47 | 153 | |
154 | // create an object that represents your choice of jet algorithm, and | |
155 | // the associated parameters | |
156 | double Rparam = header->GetRparam(); | |
157 | fastjet::Strategy strategy = header->GetStrategy(); | |
158 | fastjet::RecombinationScheme recomb_scheme = header->GetRecombScheme(); | |
159 | fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); | |
160 | fastjet::JetDefinition jet_def(algorithm, Rparam, recomb_scheme, strategy); | |
8838ab7a | 161 | |
392c9b47 | 162 | // create an object that specifies how we to define the area |
163 | fastjet::AreaDefinition area_def; | |
164 | double ghost_etamax = header->GetGhostEtaMax(); | |
165 | double ghost_area = header->GetGhostArea(); | |
166 | int active_area_repeats = header->GetActiveAreaRepeats(); | |
167 | ||
168 | // now create the object that holds info about ghosts | |
169 | fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, ghost_area); | |
170 | // and from that get an area definition | |
171 | fastjet::AreaType area_type = header->GetAreaType(); | |
172 | area_def = fastjet::AreaDefinition(area_type,ghost_spec); | |
173 | ||
8838ab7a | 174 | if(bgMode) // BG subtraction |
175 | { | |
176 | //***************************** JETS FINDING AND EXTRACTION | |
177 | // run the jet clustering with the above jet definition | |
178 | fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def); | |
392c9b47 | 179 | |
8838ab7a | 180 | // save a comment in the header |
181 | ||
182 | TString comment = "Running FastJet algorithm with the following setup. "; | |
183 | comment+= "Jet definition: "; | |
184 | comment+= TString(jet_def.description()); | |
185 | comment+= ". Area definition: "; | |
186 | comment+= TString(area_def.description()); | |
187 | comment+= ". Strategy adopted by FastJet: "; | |
188 | comment+= TString(clust_seq.strategy_string()); | |
189 | header->SetComment(comment); | |
190 | if(debug){ | |
191 | cout << "--------------------------------------------------------" << endl; | |
192 | cout << comment << endl; | |
193 | cout << "--------------------------------------------------------" << endl; | |
194 | } | |
195 | //header->PrintParameters(); | |
196 | ||
197 | ||
198 | // extract the inclusive jets with pt > ptmin, sorted by pt | |
199 | double ptmin = header->GetPtMin(); | |
200 | vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin); | |
201 | ||
202 | //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; | |
203 | ||
204 | ||
205 | //subtract background // =========================================== | |
206 | // set the rapididty , phi range within which to study the background | |
207 | double rap_max = header->GetRapMax(); | |
208 | double rap_min = header->GetRapMin(); | |
209 | double phi_max = header->GetPhiMax(); | |
210 | double phi_min = header->GetPhiMin(); | |
211 | fastjet::RangeDefinition range(rap_min, rap_max, phi_min, phi_max); | |
212 | ||
213 | // subtract background | |
214 | vector<fastjet::PseudoJet> sub_jets = clust_seq.subtracted_jets(range,ptmin); | |
215 | ||
216 | // print out | |
217 | //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n"; | |
218 | //cout << "---------------------------------------\n"; | |
219 | //cout << endl; | |
220 | //printf(" ijet rap phi Pt area +- err\n"); | |
221 | ||
222 | // sort jets into increasing pt | |
223 | vector<fastjet::PseudoJet> jets = sorted_by_pt(sub_jets); | |
224 | for (size_t j = 0; j < jets.size(); j++) { // loop for jets | |
225 | ||
226 | double area = clust_seq.area(jets[j]); | |
227 | double area_error = clust_seq.area_error(jets[j]); | |
228 | ||
229 | 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, area_error); | |
392c9b47 | 230 | |
231 | // go to write AOD info | |
8838ab7a | 232 | AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); |
233 | //cout << "Printing jet " << endl; | |
234 | if(debug) aodjet.Print(""); | |
235 | //cout << "Adding jet ... " ; | |
236 | AddJet(aodjet); | |
237 | //cout << "added \n" << endl; | |
238 | ||
239 | } | |
240 | } | |
241 | else { // No BG subtraction | |
392c9b47 | 242 | |
8838ab7a | 243 | fastjet::ClusterSequence clust_seq(input_particles, jet_def); |
244 | ||
245 | // save a comment in the header | |
246 | ||
247 | TString comment = "Running FastJet algorithm with the following setup. "; | |
248 | comment+= "Jet definition: "; | |
249 | comment+= TString(jet_def.description()); | |
250 | comment+= ". Strategy adopted by FastJet: "; | |
251 | comment+= TString(clust_seq.strategy_string()); | |
252 | header->SetComment(comment); | |
253 | if(debug){ | |
254 | cout << "--------------------------------------------------------" << endl; | |
255 | cout << comment << endl; | |
256 | cout << "--------------------------------------------------------" << endl; | |
257 | } | |
258 | //header->PrintParameters(); | |
259 | ||
260 | // extract the inclusive jets with pt > ptmin, sorted by pt | |
261 | double ptmin = header->GetPtMin(); | |
262 | vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin); | |
263 | ||
264 | //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; | |
392c9b47 | 265 | |
8838ab7a | 266 | vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusive_jets); // Added by me |
267 | for (size_t j = 0; j < jets.size(); j++) { // loop for jets // Added by me | |
268 | ||
269 | printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp()); | |
270 | ||
271 | // go to write AOD info | |
272 | AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); | |
273 | //cout << "Printing jet " << endl; | |
274 | if(debug) aodjet.Print(""); | |
275 | //cout << "Adding jet ... " ; | |
276 | AddJet(aodjet); | |
277 | //cout << "added \n" << endl; | |
278 | ||
279 | } // end loop for jets | |
280 | } | |
281 | ||
a17e6965 | 282 | } |
283 | ||
392c9b47 | 284 | //____________________________________________________________________________ |
285 | void AliFastJetFinder::RunTest(const char* datafile) | |
a17e6965 | 286 | |
a17e6965 | 287 | { |
a17e6965 | 288 | |
392c9b47 | 289 | // This simple test run the kt algorithm for an ascii file testdata.dat |
290 | // read input particles ----------------------------- | |
291 | vector<fastjet::PseudoJet> input_particles; | |
292 | Float_t px,py,pz,en; | |
293 | ifstream in; | |
294 | Int_t nlines = 0; | |
295 | // we assume a file basic.dat in the current directory | |
296 | // this file has 3 columns of float data | |
297 | in.open(datafile); | |
298 | while (1) { | |
299 | in >> px >> py >> pz >> en; | |
300 | if (!in.good()) break; | |
301 | //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz); | |
302 | nlines++; | |
303 | input_particles.push_back(fastjet::PseudoJet(px,py,pz,en)); | |
a17e6965 | 304 | } |
392c9b47 | 305 | //printf(" found %d pointsn",nlines); |
306 | in.close(); | |
307 | ////////////////////////////////////////////////// | |
308 | ||
309 | // create an object that represents your choice of jet algorithm, and | |
310 | // the associated parameters | |
311 | double Rparam = 1.0; | |
312 | fastjet::Strategy strategy = fastjet::Best; | |
313 | fastjet::RecombinationScheme recomb_scheme = fastjet::BIpt_scheme; | |
314 | fastjet::JetDefinition jet_def(fastjet::kt_algorithm, Rparam, recomb_scheme, strategy); | |
315 | ||
316 | ||
317 | ||
318 | // create an object that specifies how we to define the area | |
319 | fastjet::AreaDefinition area_def; | |
320 | double ghost_etamax = 7.0; | |
321 | double ghost_area = 0.05; | |
322 | int active_area_repeats = 1; | |
323 | ||
324 | ||
325 | // now create the object that holds info about ghosts | |
326 | fastjet::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, ghost_area); | |
327 | // and from that get an area definition | |
328 | area_def = fastjet::AreaDefinition(fastjet::active_area,ghost_spec); | |
329 | ||
330 | ||
331 | // run the jet clustering with the above jet definition | |
332 | fastjet::ClusterSequenceArea clust_seq(input_particles, jet_def, area_def); | |
333 | ||
334 | ||
335 | // tell the user what was done | |
336 | cout << "--------------------------------------------------------" << endl; | |
337 | cout << "Jet definition was: " << jet_def.description() << endl; | |
338 | cout << "Area definition was: " << area_def.description() << endl; | |
339 | cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl; | |
340 | cout << "--------------------------------------------------------" << endl; | |
341 | ||
342 | ||
343 | // extract the inclusive jets with pt > 5 GeV, sorted by pt | |
344 | double ptmin = 5.0; | |
345 | vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin); | |
346 | ||
347 | cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; | |
348 | ||
349 | ||
350 | //subtract background // =========================================== | |
351 | // set the rapididty range within which to study the background | |
352 | double rap_max = ghost_etamax - Rparam; | |
353 | fastjet::RangeDefinition range(rap_max); | |
354 | // subtract background | |
355 | vector<fastjet::PseudoJet> sub_jets = clust_seq.subtracted_jets(range,ptmin); | |
356 | ||
357 | // print them out //================================================ | |
358 | cout << "Printing inclusive jets after background subtraction \n"; | |
359 | cout << "------------------------------------------------------\n"; | |
360 | // sort jets into increasing pt | |
361 | vector<fastjet::PseudoJet> jets = sorted_by_pt(sub_jets); | |
362 | ||
363 | printf(" ijet rap phi Pt area +- err\n"); | |
364 | for (size_t j = 0; j < jets.size(); j++) { | |
365 | ||
366 | double area = clust_seq.area(jets[j]); | |
367 | double area_error = clust_seq.area_error(jets[j]); | |
368 | ||
8838ab7a | 369 | printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(), |
392c9b47 | 370 | jets[j].phi(),jets[j].perp(), area, area_error); |
371 | } | |
372 | cout << endl; | |
373 | // ================================================================ | |
a17e6965 | 374 | |
392c9b47 | 375 | |
376 | ||
a17e6965 | 377 | } |
378 | ||
392c9b47 | 379 | //____________________________________________________________________________ |
a17e6965 | 380 | |
381 | void AliFastJetFinder::WriteJHeaderToFile() | |
382 | { | |
a17e6965 | 383 | fHeader->Write(); |
384 | } | |
385 | ||
9157b9c9 | 386 | //____________________________________________________________________________ |
8838ab7a | 387 | |
388 | Float_t AliFastJetFinder::EtaToTheta(Float_t arg) | |
389 | { | |
390 | // return (180./TMath::Pi())*2.*atan(exp(-arg)); | |
391 | return 2.*atan(exp(-arg)); | |
392 | ||
393 | ||
394 | } | |
395 | ||
396 | //____________________________________________________________________________ | |
397 | ||
398 | void AliFastJetFinder::InitTask(TChain *tree) | |
399 | { | |
400 | ||
401 | printf("Fast jet finder initialization ******************"); | |
402 | fReader->CreateTasks(tree); | |
403 | ||
404 | } |