1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //---------------------------------------------------------------------
19 // FastJet v2.3.4 finder algorithm interface
20 // Last modification: Neutral cell energy included in the jet reconstruction
22 // Authors: Rafael.Diaz.Valdes@cern.ch
23 // Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option)
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
28 //---------------------------------------------------------------------
30 #include <Riostream.h>
32 #include "AliFastJetFinder.h"
33 #include "AliFastJetHeaderV1.h"
34 #include "AliFastJetInput.h"
35 #include "AliFastJetBkg.h"
36 #include "AliAODJetEventBackground.h"
37 #include "AliAODJet.h"
39 #include "fastjet/PseudoJet.hh"
40 #include "fastjet/ClusterSequenceArea.hh"
41 #include "fastjet/AreaDefinition.hh"
42 #include "fastjet/JetDefinition.hh"
48 ClassImp(AliFastJetFinder)
50 ////////////////////////////////////////////////////////////////////////
52 AliFastJetFinder::AliFastJetFinder():
54 fInputFJ(new AliFastJetInput()),
55 fJetBkg(new AliFastJetBkg())
60 //____________________________________________________________________________
61 AliFastJetFinder::~AliFastJetFinder()
69 //______________________________________________________________________________
70 void AliFastJetFinder::FindJets()
72 // runs a FASTJET based algo
74 //pick up fastjet header
75 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
76 Int_t debug = header->GetDebug(); // debug option
77 Bool_t bgMode = header->GetBGMode(); // choose to subtract BG or not
78 if(debug>0) cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl;
81 // read input particles -----------------------------
83 vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
84 if(inputParticles.size()==0){
85 if(debug>0) Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
89 // create an object that represents your choice of jet algorithm, and
90 // the associated parameters
91 double rParam = header->GetRparam();
92 double rBkgParam = header->GetRparamBkg();
93 fastjet::Strategy strategy = header->GetStrategy();
94 fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
95 fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
96 fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy);
98 // create an object that specifies how we to define the area
99 fastjet::AreaDefinition areaDef;
100 double ghostEtamax = header->GetGhostEtaMax();
101 double ghostArea = header->GetGhostArea();
102 int activeAreaRepeats = header->GetActiveAreaRepeats();
104 // now create the object that holds info about ghosts
105 fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
106 // and from that get an area definition
107 fastjet::AreaType areaType = header->GetAreaType();
108 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
110 //***************************** JETS FINDING
111 // run the jet clustering with the above jet definition
112 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
114 vector<fastjet::PseudoJet> jets;
116 if(bgMode) // Do BG subtraction directly with the same algorithm (cambridge or kt) for jet signal and background
118 //***************************** JETS FINDING FOR RHO ESTIMATION
119 // run the jet clustering with the above jet definition
120 fastjet::JetAlgorithm algorithmBkg = header->GetBGAlgorithm();
121 fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy);
122 fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef);
124 // save a comment in the header
125 TString comment = "Running FastJet algorithm with the following setup. ";
126 comment+= "Jet definition: ";
127 comment+= TString(jetDef.description());
128 comment+= "Jet bckg definition: ";
129 comment+= TString(jetDefBkg.description());
130 comment+= ". Area definition: ";
131 comment+= TString(areaDef.description());
132 comment+= ". Strategy adopted by FastJet: ";
133 comment+= TString(clust_seq.strategy_string());
134 header->SetComment(comment);
136 cout << "--------------------------------------------------------" << endl;
137 cout << comment << endl;
138 cout << "--------------------------------------------------------" << endl;
141 // extract the inclusive jets sorted by pt
142 double ptmin = header->GetPtMin();
143 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets();
145 //subtract background // ===========================================
146 // set the rapididty , phi range within which to study the background
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);
153 // Extract rho and sigma
156 Double_t meanarea = 0.;
157 Bool_t kUse4VectorArea = header->Use4VectorArea();
158 vector<fastjet::PseudoJet> bkgJets = clust_seq_bkg.inclusive_jets();
159 clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, kUse4VectorArea, rho, sigma, meanarea, false);
161 // subtract background and extract jets bkg subtracted
162 vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptmin);
164 // sort jets into increasing pt
165 jets = sorted_by_pt(subJets);
169 else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
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);
179 cout << "--------------------------------------------------------" << endl;
180 cout << comment << endl;
181 cout << "--------------------------------------------------------" << endl;
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);
188 jets = sorted_by_pt(inclusiveJets);
192 for (size_t j = 0; j < jets.size(); j++) { // loop for jets
194 double area = clust_seq.area(jets[j]);
195 double areaError = clust_seq.area_error(jets[j]);
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());
199 vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
200 int nCon= constituents.size();
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
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("");
215 for (int i=0; i < nCon; i++)
217 fastjet::PseudoJet mPart=constituents[i];
218 ind[i]=mPart.user_index();
220 // Jet constituents (charged tracks) added to the AliAODJet
221 AliJetCalTrkEvent* calEvt = GetCalTrkEvent();
222 for(Int_t itrack=0; itrack<calEvt->GetNCalTrkTracks(); itrack++)
226 TObject *track = calEvt->GetCalTrkTrack(itrack)->GetTrackObject();
227 aodjet.AddTrack(track);
230 } // End loop on Constituents
234 } // end loop for jets
238 //____________________________________________________________________________
239 void AliFastJetFinder::RunTest(const char* datafile)
241 // This simple test run the kt algorithm for an ascii file testdata.dat
243 // read input particles -----------------------------
244 vector<fastjet::PseudoJet> inputParticles;
248 // we assume a file basic.dat in the current directory
249 // this file has 3 columns of float data
252 in >> px >> py >> pz >> en;
253 if (!in.good()) break;
254 //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz);
256 inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en));
258 //printf(" found %d pointsn",nlines);
260 //////////////////////////////////////////////////
262 // create an object that represents your choice of jet algorithm, and
263 // the associated parameters
265 fastjet::Strategy strategy = fastjet::Best;
266 fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme;
267 fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy);
269 // create an object that specifies how we to define the area
270 fastjet::AreaDefinition areaDef;
271 double ghostEtamax = 7.0;
272 double ghostArea = 0.05;
273 int activeAreaRepeats = 1;
275 // now create the object that holds info about ghosts
276 fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
277 // and from that get an area definition
278 areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec);
280 // run the jet clustering with the above jet definition
281 fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
283 // tell the user what was done
284 cout << "--------------------------------------------------------" << endl;
285 cout << "Jet definition was: " << jetDef.description() << endl;
286 cout << "Area definition was: " << areaDef.description() << endl;
287 cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl;
288 cout << "--------------------------------------------------------" << endl;
290 // extract the inclusive jets with pt > 5 GeV, sorted by pt
292 vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
294 cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
296 //subtract background // ===========================================
297 // set the rapididty range within which to study the background
298 double rapMax = ghostEtamax - rParam;
299 fastjet::RangeDefinition range(rapMax);
300 // subtract background
301 vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(range,ptmin);
303 // print them out //================================================
304 cout << "Printing inclusive jets after background subtraction \n";
305 cout << "------------------------------------------------------\n";
306 // sort jets into increasing pt
307 vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
309 printf(" ijet rap phi Pt area +- err\n");
310 for (size_t j = 0; j < jets.size(); j++) {
311 double area = clust_seq.area(jets[j]);
312 double areaError = clust_seq.area_error(jets[j]);
313 printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(),
314 jets[j].phi(),jets[j].perp(), area, areaError);
317 // ================================================================
321 //____________________________________________________________________________
322 void AliFastJetFinder::WriteJHeaderToFile() const
324 // Write Jet Header to file
329 //____________________________________________________________________________
330 Bool_t AliFastJetFinder::ProcessEvent()
333 // Charged only or charged+neutral jets
335 fInputFJ->SetHeader(fHeader);
336 fInputFJ->SetCalTrkEvent(GetCalTrkEvent());
337 fInputFJ->FillInput();
344 AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
346 fJetBkg->SetHeader(fHeader);
347 fJetBkg->SetFastJetInput(fInputFJ);
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);
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);