]>
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 | ||
139cbd96 | 16 | /* $Id$ */ |
392c9b47 | 17 | |
a17e6965 | 18 | //--------------------------------------------------------------------- |
392c9b47 | 19 | // FastJet v2.3.4 finder algorithm interface |
8838ab7a | 20 | // Last modification: Neutral cell energy included in the jet reconstruction |
21 | // | |
22 | // Authors: Rafael.Diaz.Valdes@cern.ch | |
23 | // Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option) | |
392c9b47 | 24 | // |
139cbd96 | 25 | // ** 2011 magali.estienne@subatech.in2p3.fr & alexandre.shabetai@cern.ch |
26 | // new implementation of background subtraction | |
27 | // allowing to subtract bkg using a different algo than the one used for signal jets | |
a17e6965 | 28 | //--------------------------------------------------------------------- |
29 | ||
30 | #include <Riostream.h> | |
a17e6965 | 31 | |
32 | #include "AliFastJetFinder.h" | |
8838ab7a | 33 | #include "AliFastJetHeaderV1.h" |
d1993270 | 34 | #include "AliFastJetInput.h" |
139cbd96 | 35 | #include "AliFastJetBkg.h" |
b8bf1e90 | 36 | #include "AliAODJetEventBackground.h" |
139cbd96 | 37 | #include "AliAODJet.h" |
392c9b47 | 38 | |
39 | #include "fastjet/PseudoJet.hh" | |
40 | #include "fastjet/ClusterSequenceArea.hh" | |
41 | #include "fastjet/AreaDefinition.hh" | |
42 | #include "fastjet/JetDefinition.hh" | |
392c9b47 | 43 | |
392c9b47 | 44 | #include<vector> |
392c9b47 | 45 | |
46 | using namespace std; | |
a17e6965 | 47 | |
9157b9c9 | 48 | ClassImp(AliFastJetFinder) |
8838ab7a | 49 | |
139cbd96 | 50 | //////////////////////////////////////////////////////////////////////// |
a17e6965 | 51 | |
52 | AliFastJetFinder::AliFastJetFinder(): | |
856618e7 | 53 | AliJetFinder(), |
287697fc | 54 | fInputFJ(new AliFastJetInput()), |
139cbd96 | 55 | fJetBkg(new AliFastJetBkg()) |
a17e6965 | 56 | { |
57 | // Constructor | |
58 | } | |
59 | ||
392c9b47 | 60 | //____________________________________________________________________________ |
a17e6965 | 61 | AliFastJetFinder::~AliFastJetFinder() |
a17e6965 | 62 | { |
63 | // destructor | |
139cbd96 | 64 | delete fInputFJ; |
65 | delete fJetBkg; | |
66 | ||
a17e6965 | 67 | } |
68 | ||
392c9b47 | 69 | //______________________________________________________________________________ |
a17e6965 | 70 | void AliFastJetFinder::FindJets() |
a17e6965 | 71 | { |
139cbd96 | 72 | // runs a FASTJET based algo |
d56792bf | 73 | |
392c9b47 | 74 | //pick up fastjet header |
8838ab7a | 75 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; |
cc6a2227 | 76 | Int_t debug = header->GetDebug(); // debug option |
8838ab7a | 77 | Bool_t bgMode = header->GetBGMode(); // choose to subtract BG or not |
139cbd96 | 78 | if(debug>0) cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl; |
8f1fc92d | 79 | |
392c9b47 | 80 | // RUN ALGORITHM |
81 | // read input particles ----------------------------- | |
d1993270 | 82 | |
83 | vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles(); | |
03ecc336 | 84 | if(inputParticles.size()==0){ |
139cbd96 | 85 | if(debug>0) Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__); |
03ecc336 | 86 | return; |
87 | } | |
d1993270 | 88 | |
392c9b47 | 89 | // create an object that represents your choice of jet algorithm, and |
90 | // the associated parameters | |
3dda898d | 91 | double rParam = header->GetRparam(); |
8f1fc92d | 92 | double rBkgParam = header->GetRparamBkg(); |
392c9b47 | 93 | fastjet::Strategy strategy = header->GetStrategy(); |
3dda898d | 94 | fastjet::RecombinationScheme recombScheme = header->GetRecombScheme(); |
392c9b47 | 95 | fastjet::JetAlgorithm algorithm = header->GetAlgorithm(); |
d1993270 | 96 | fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy); |
8838ab7a | 97 | |
392c9b47 | 98 | // create an object that specifies how we to define the area |
3dda898d | 99 | fastjet::AreaDefinition areaDef; |
139cbd96 | 100 | double ghostEtamax = header->GetGhostEtaMax(); |
101 | double ghostArea = header->GetGhostArea(); | |
3dda898d | 102 | int activeAreaRepeats = header->GetActiveAreaRepeats(); |
392c9b47 | 103 | |
104 | // now create the object that holds info about ghosts | |
d1993270 | 105 | fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea); |
392c9b47 | 106 | // and from that get an area definition |
3dda898d | 107 | fastjet::AreaType areaType = header->GetAreaType(); |
d1993270 | 108 | areaDef = fastjet::AreaDefinition(areaType,ghostSpec); |
392c9b47 | 109 | |
8f1fc92d | 110 | //***************************** JETS FINDING |
111 | // run the jet clustering with the above jet definition | |
112 | fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef); | |
a38beabd | 113 | |
8f1fc92d | 114 | vector<fastjet::PseudoJet> jets; |
a38beabd | 115 | |
8f1fc92d | 116 | if(bgMode) // Do BG subtraction directly with the same algorithm (cambridge or kt) for jet signal and background |
8838ab7a | 117 | { |
139cbd96 | 118 | //***************************** JETS FINDING FOR RHO ESTIMATION |
8838ab7a | 119 | // run the jet clustering with the above jet definition |
139cbd96 | 120 | fastjet::JetAlgorithm algorithmBkg = header->GetBGAlgorithm(); |
121 | fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy); | |
8f1fc92d | 122 | fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef); |
392c9b47 | 123 | |
8838ab7a | 124 | // save a comment in the header |
8838ab7a | 125 | TString comment = "Running FastJet algorithm with the following setup. "; |
126 | comment+= "Jet definition: "; | |
d1993270 | 127 | comment+= TString(jetDef.description()); |
8f1fc92d | 128 | comment+= "Jet bckg definition: "; |
129 | comment+= TString(jetDefBkg.description()); | |
8838ab7a | 130 | comment+= ". Area definition: "; |
3dda898d | 131 | comment+= TString(areaDef.description()); |
139cbd96 | 132 | comment+= ". Strategy adopted by FastJet: "; |
8838ab7a | 133 | comment+= TString(clust_seq.strategy_string()); |
134 | header->SetComment(comment); | |
139cbd96 | 135 | if(debug>0){ |
136 | cout << "--------------------------------------------------------" << endl; | |
137 | cout << comment << endl; | |
138 | cout << "--------------------------------------------------------" << endl; | |
8838ab7a | 139 | } |
8838ab7a | 140 | |
139cbd96 | 141 | // extract the inclusive jets sorted by pt |
8838ab7a | 142 | double ptmin = header->GetPtMin(); |
8f1fc92d | 143 | vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(); |
139cbd96 | 144 | |
8838ab7a | 145 | //subtract background // =========================================== |
146 | // set the rapididty , phi range within which to study the background | |
3dda898d | 147 | double rapMax = header->GetRapMax(); |
148 | double rapMin = header->GetRapMin(); | |
149 | double phiMax = header->GetPhiMax(); | |
150 | double phiMin = header->GetPhiMin(); | |
151 | fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax); | |
139cbd96 | 152 | |
8f1fc92d | 153 | // Extract rho and sigma |
154 | Double_t rho = 0.; | |
155 | Double_t sigma = 0.; | |
156 | Double_t meanarea = 0.; | |
139cbd96 | 157 | Bool_t kUse4VectorArea = header->Use4VectorArea(); |
8f1fc92d | 158 | vector<fastjet::PseudoJet> bkgJets = clust_seq_bkg.inclusive_jets(); |
139cbd96 | 159 | clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, kUse4VectorArea, rho, sigma, meanarea, false); |
8f1fc92d | 160 | |
139cbd96 | 161 | // subtract background and extract jets bkg subtracted |
8f1fc92d | 162 | vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptmin); |
8838ab7a | 163 | |
164 | // sort jets into increasing pt | |
8f1fc92d | 165 | jets = sorted_by_pt(subJets); |
166 | ||
8838ab7a | 167 | } |
a38beabd | 168 | |
139cbd96 | 169 | else { // No BG subtraction!!!!!!!! Default header is bgmode=0. |
170 | ||
171 | // save a comment in the header | |
172 | TString comment = "Running FastJet algorithm with the following setup. "; | |
173 | comment+= "Jet definition: "; | |
174 | comment+= TString(jetDef.description()); | |
175 | comment+= ". Strategy adopted by FastJet: "; | |
176 | comment+= TString(clust_seq.strategy_string()); | |
177 | header->SetComment(comment); | |
178 | if(debug>0){ | |
179 | cout << "--------------------------------------------------------" << endl; | |
180 | cout << comment << endl; | |
181 | cout << "--------------------------------------------------------" << endl; | |
182 | } | |
183 | ||
184 | // extract the inclusive jets with pt > ptmin, sorted by pt | |
185 | double ptmin = header->GetPtMin(); | |
186 | vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin); | |
187 | ||
188 | jets = sorted_by_pt(inclusiveJets); | |
189 | ||
8f1fc92d | 190 | } |
139cbd96 | 191 | |
192 | for (size_t j = 0; j < jets.size(); j++) { // loop for jets | |
a38beabd | 193 | |
139cbd96 | 194 | double area = clust_seq.area(jets[j]); |
195 | double areaError = clust_seq.area_error(jets[j]); | |
392c9b47 | 196 | |
139cbd96 | 197 | if(debug>0) printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp()); |
d1993270 | 198 | |
139cbd96 | 199 | vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]); |
200 | int nCon= constituents.size(); | |
201 | TArrayI ind(nCon); | |
202 | ||
203 | if ((jets[j].eta() > (header->GetJetEtaMax())) || | |
204 | (jets[j].eta() < (header->GetJetEtaMin())) || | |
205 | (jets[j].phi() > (header->GetJetPhiMax())) || | |
206 | (jets[j].phi() < (header->GetJetPhiMin())) || | |
207 | (jets[j].perp() < header->GetPtMin())) continue; // acceptance eta range and etmin | |
208 | ||
209 | // go to write AOD info | |
210 | AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E()); | |
211 | aodjet.SetEffArea(area,areaError); | |
212 | //cout << "Printing jet " << endl; | |
213 | if(debug>0) aodjet.Print(""); | |
8f1fc92d | 214 | |
139cbd96 | 215 | for (int i=0; i < nCon; i++) |
216 | { | |
217 | fastjet::PseudoJet mPart=constituents[i]; | |
218 | ind[i]=mPart.user_index(); | |
219 | ||
220 | // Jet constituents (charged tracks) added to the AliAODJet | |
221 | AliJetCalTrkEvent* calEvt = GetCalTrkEvent(); | |
222 | for(Int_t itrack=0; itrack<calEvt->GetNCalTrkTracks(); itrack++) | |
d1993270 | 223 | { |
139cbd96 | 224 | if(itrack==ind[i]) |
225 | { | |
226 | TObject *track = calEvt->GetCalTrkTrack(itrack)->GetTrackObject(); | |
227 | aodjet.AddTrack(track); | |
228 | } | |
229 | } | |
230 | } // End loop on Constituents | |
231 | ||
232 | AddJet(aodjet); | |
233 | ||
234 | } // end loop for jets | |
8838ab7a | 235 | |
a17e6965 | 236 | } |
237 | ||
392c9b47 | 238 | //____________________________________________________________________________ |
239 | void AliFastJetFinder::RunTest(const char* datafile) | |
a17e6965 | 240 | { |
139cbd96 | 241 | // This simple test run the kt algorithm for an ascii file testdata.dat |
a17e6965 | 242 | |
139cbd96 | 243 | // read input particles ----------------------------- |
3dda898d | 244 | vector<fastjet::PseudoJet> inputParticles; |
392c9b47 | 245 | Float_t px,py,pz,en; |
246 | ifstream in; | |
247 | Int_t nlines = 0; | |
248 | // we assume a file basic.dat in the current directory | |
249 | // this file has 3 columns of float data | |
250 | in.open(datafile); | |
251 | while (1) { | |
139cbd96 | 252 | in >> px >> py >> pz >> en; |
253 | if (!in.good()) break; | |
254 | //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz); | |
255 | nlines++; | |
256 | inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en)); | |
257 | } | |
258 | //printf(" found %d pointsn",nlines); | |
259 | in.close(); | |
260 | ////////////////////////////////////////////////// | |
392c9b47 | 261 | |
262 | // create an object that represents your choice of jet algorithm, and | |
263 | // the associated parameters | |
3dda898d | 264 | double rParam = 1.0; |
392c9b47 | 265 | fastjet::Strategy strategy = fastjet::Best; |
3dda898d | 266 | fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme; |
d1993270 | 267 | fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy); |
392c9b47 | 268 | |
392c9b47 | 269 | // create an object that specifies how we to define the area |
3dda898d | 270 | fastjet::AreaDefinition areaDef; |
271 | double ghostEtamax = 7.0; | |
272 | double ghostArea = 0.05; | |
273 | int activeAreaRepeats = 1; | |
392c9b47 | 274 | |
392c9b47 | 275 | // now create the object that holds info about ghosts |
d1993270 | 276 | fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea); |
392c9b47 | 277 | // and from that get an area definition |
d1993270 | 278 | areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec); |
392c9b47 | 279 | |
280 | // run the jet clustering with the above jet definition | |
d1993270 | 281 | fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef); |
392c9b47 | 282 | |
392c9b47 | 283 | // tell the user what was done |
284 | cout << "--------------------------------------------------------" << endl; | |
d1993270 | 285 | cout << "Jet definition was: " << jetDef.description() << endl; |
3dda898d | 286 | cout << "Area definition was: " << areaDef.description() << endl; |
392c9b47 | 287 | cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl; |
288 | cout << "--------------------------------------------------------" << endl; | |
392c9b47 | 289 | |
290 | // extract the inclusive jets with pt > 5 GeV, sorted by pt | |
291 | double ptmin = 5.0; | |
3dda898d | 292 | vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin); |
392c9b47 | 293 | |
294 | cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl; | |
295 | ||
392c9b47 | 296 | //subtract background // =========================================== |
297 | // set the rapididty range within which to study the background | |
3dda898d | 298 | double rapMax = ghostEtamax - rParam; |
299 | fastjet::RangeDefinition range(rapMax); | |
392c9b47 | 300 | // subtract background |
3dda898d | 301 | vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(range,ptmin); |
392c9b47 | 302 | |
303 | // print them out //================================================ | |
304 | cout << "Printing inclusive jets after background subtraction \n"; | |
305 | cout << "------------------------------------------------------\n"; | |
306 | // sort jets into increasing pt | |
3dda898d | 307 | vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets); |
392c9b47 | 308 | |
309 | printf(" ijet rap phi Pt area +- err\n"); | |
310 | for (size_t j = 0; j < jets.size(); j++) { | |
139cbd96 | 311 | double area = clust_seq.area(jets[j]); |
3dda898d | 312 | double areaError = clust_seq.area_error(jets[j]); |
8838ab7a | 313 | printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(), |
3dda898d | 314 | jets[j].phi(),jets[j].perp(), area, areaError); |
392c9b47 | 315 | } |
316 | cout << endl; | |
317 | // ================================================================ | |
a17e6965 | 318 | |
a17e6965 | 319 | } |
320 | ||
392c9b47 | 321 | //____________________________________________________________________________ |
b430dff0 | 322 | void AliFastJetFinder::WriteJHeaderToFile() const |
a17e6965 | 323 | { |
139cbd96 | 324 | // Write Jet Header to file |
a17e6965 | 325 | fHeader->Write(); |
8838ab7a | 326 | |
327 | } | |
328 | ||
329 | //____________________________________________________________________________ | |
287697fc | 330 | Bool_t AliFastJetFinder::ProcessEvent() |
331 | { | |
287697fc | 332 | // Process one event |
139cbd96 | 333 | // Charged only or charged+neutral jets |
287697fc | 334 | |
287697fc | 335 | fInputFJ->SetHeader(fHeader); |
139cbd96 | 336 | fInputFJ->SetCalTrkEvent(GetCalTrkEvent()); |
287697fc | 337 | fInputFJ->FillInput(); |
139cbd96 | 338 | |
339 | // Find Jets | |
287697fc | 340 | FindJets(); |
341 | ||
139cbd96 | 342 | // Background |
343 | if(fAODEvBkg){ | |
344 | AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader; | |
294056ab | 345 | |
21ead515 | 346 | fJetBkg->SetHeader(fHeader); |
21ead515 | 347 | fJetBkg->SetFastJetInput(fInputFJ); |
139cbd96 | 348 | |
349 | Int_t count = 0; | |
350 | if(header->GetBkgFastJetb()){ | |
351 | Double_t sigma1 = 0 , meanarea1= 0, bkg1 = 0; | |
352 | fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1); | |
353 | fAODEvBkg->SetBackground(count,bkg1,sigma1,meanarea1); | |
354 | count++; | |
355 | } | |
356 | ||
357 | if(header->GetBkgFastJetWoHardest()){ | |
358 | Double_t sigma2 = 0 , meanarea2 = 0, bkg2 = 0; | |
359 | fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2); | |
360 | fAODEvBkg->SetBackground(count,bkg2,sigma2,meanarea2); | |
361 | count++; | |
362 | } | |
363 | ||
364 | } | |
294056ab | 365 | |
287697fc | 366 | Reset(); |
367 | return kTRUE; | |
368 | ||
369 | } |