Ensure vacuum inside the beam pipe for upgrade (Mario)
[u/mrichter/AliRoot.git] / JETAN / AliFastJetFinder.cxx
CommitLineData
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
46using namespace std;
a17e6965 47
9157b9c9 48ClassImp(AliFastJetFinder)
8838ab7a 49
139cbd96 50////////////////////////////////////////////////////////////////////////
a17e6965 51
52AliFastJetFinder::AliFastJetFinder():
856618e7 53 AliJetFinder(),
287697fc 54 fInputFJ(new AliFastJetInput()),
139cbd96 55 fJetBkg(new AliFastJetBkg())
a17e6965 56{
57 // Constructor
58}
59
392c9b47 60//____________________________________________________________________________
a17e6965 61AliFastJetFinder::~AliFastJetFinder()
a17e6965 62{
63 // destructor
139cbd96 64 delete fInputFJ;
65 delete fJetBkg;
66
a17e6965 67}
68
392c9b47 69//______________________________________________________________________________
a17e6965 70void 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//____________________________________________________________________________
239void 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 322void AliFastJetFinder::WriteJHeaderToFile() const
a17e6965 323{
139cbd96 324 // Write Jet Header to file
a17e6965 325 fHeader->Write();
8838ab7a 326
327}
328
329//____________________________________________________________________________
287697fc 330Bool_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}