Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / JETAN / AliSISConeJetFinder.cxx
CommitLineData
2551af9d 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
17//---------------------------------------------------------------------
8838ab7a 18// FastJet v2.3.4 finder algorithm interface
2551af9d 19//
20// Author: swensy.jangal@ires.in2p3.fr
8838ab7a 21//
22// Last modification: Neutral cell energy included in the jet reconstruction
23//
24// Author: Magali.estienne@subatech.in2p3.fr
2551af9d 25//---------------------------------------------------------------------
26
8838ab7a 27
2551af9d 28#include <Riostream.h>
29#include <TArrayF.h>
2551af9d 30#include <TFile.h>
31#include <TH1F.h>
32#include <TH2F.h>
33#include <TLorentzVector.h>
34#include <TRandom.h>
8838ab7a 35#include <TClonesArray.h>
2551af9d 36
37#include "AliHeader.h"
2551af9d 38#include "AliJetKineReader.h"
39#include "AliJetReader.h"
40#include "AliJetReaderHeader.h"
8838ab7a 41#include "AliJetUnitArray.h"
2551af9d 42#include "AliSISConeJetFinder.h"
43#include "AliSISConeJetHeader.h"
44
45#include "fastjet/AreaDefinition.hh"
46#include "fastjet/ClusterSequenceArea.hh"
47#include "fastjet/JetDefinition.hh"
48#include "fastjet/PseudoJet.hh"
49
8838ab7a 50// get info on how fastjet was configured
2551af9d 51#include "fastjet/config.h"
52
53#ifdef ENABLE_PLUGIN_SISCONE
54#include "fastjet/SISConePlugin.hh"
55#endif
56
8838ab7a 57#include<sstream> // needed for internal io
2551af9d 58#include<vector>
59#include <cmath>
60
61using namespace std;
62
8838ab7a 63
2551af9d 64ClassImp(AliSISConeJetFinder)
65
8838ab7a 66
2551af9d 67//____________________________________________________________________________
68
69AliSISConeJetFinder::AliSISConeJetFinder():
70AliJetFinder()
71{
72 // Constructor
73}
74
75//____________________________________________________________________________
76
77AliSISConeJetFinder::~AliSISConeJetFinder()
78{
79 // destructor
80}
81
82//______________________________________________________________________________
83void AliSISConeJetFinder::FindJets()
84{
85
2551af9d 86 // Pick up siscone header
73faae2f 87 AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader;
88 AliJetReaderHeader *readerheader = (AliJetReaderHeader*)fReader->GetReaderHeader();
89
90 Int_t debug = header->GetDebug(); // debug option
91 Int_t fOpt = fReader->GetReaderHeader()->GetDetector();
2551af9d 92
2551af9d 93
94//******************************** SISCONE PLUGIN CONFIGURATION
95// Here we look for SISCone parameters in the header and we define our plugin.
96
73faae2f 97 Double_t coneRadius = header->GetConeRadius(); // cone radius
98 Double_t overlapThreshold = header->GetOverlapThreshold(); // overlap parameter
99 Int_t nPassMax = header->GetNPassMax(); // maximum number of passes
100 Double_t ptProtoJetMin = header->GetPtProtojetMin(); // pT min of protojets
101 Double_t caching = header->GetCaching(); // do we record found cones for this set of data?
2551af9d 102
8838ab7a 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
2551af9d 105
106 fastjet::JetDefinition::Plugin * plugin;
107 plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
108
8838ab7a 109 vector<fastjet::PseudoJet> inputParticles;
73faae2f 110
2551af9d 111//******************************** READING OF INPUT PARTICLES
112// Here we look for px, py pz and energy of each particle that we gather in a PseudoJet object, and we put all these PseudoJet in a vector of PseudoJets : input_particles.
113
73faae2f 114 Double_t pttrackcut = readerheader->GetPtCut();
115 if (debug)cout<<"pT track cut for SISCone jet finder = "<<pttrackcut<<" GeV"<<endl;
489f77cd 116
73faae2f 117 if(fOpt==0)
118 {
119 TClonesArray *lvArray = fReader->GetMomentumArray();
120
121 // We check if lvArray is ok
122 if(lvArray == 0)
123 {
124 cout << "Could not get the momentum array" << endl;
125 delete plugin;
126 return;
127 }
8838ab7a 128
73faae2f 129 Int_t nIn = lvArray->GetEntries();
489f77cd 130
73faae2f 131 if(nIn == 0)// nIn = Number of particles in the event
132 {
133 if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
134 delete plugin;
135 return;
136 }
8838ab7a 137
73faae2f 138 Float_t px,py,pz,en;
8838ab7a 139
73faae2f 140 // Load input vectors
141 for(Int_t i = 0; i < nIn; i++)
142 {
143 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
144 px = lv->Px();
145 py = lv->Py();
146 pz = lv->Pz();
147 en = lv->Energy();
148
149 if (pttrackcut > TMath::Sqrt(px*px+py*py)) continue; //added by syssy
8838ab7a 150
73faae2f 151 fastjet::PseudoJet inputPart(px,py,pz,en);
152 inputPart.set_user_index(i);
153 inputParticles.push_back(inputPart);
154 }
155 }
156 else
157 {
158 TClonesArray* fUnit = fReader->GetUnitArray();
159 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; delete plugin; return; }
160 Int_t nIn = fUnit->GetEntries();
161 if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; delete plugin; return; }
162 // Information extracted from fUnitArray
163 // load input vectors and calculate total energy in array
164 Float_t pt,eta,phi,theta,px,py,pz,en;
165 Int_t ipart = 0;
166 for(Int_t i=0; i<nIn; i++)
167 {
168 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
8838ab7a 169
73faae2f 170 if(uArray->GetUnitEnergy()>0.)
171 {
172 // It is not necessary anymore to cut on particle pt
173 pt = uArray->GetUnitEnergy();
174 eta = uArray->GetUnitEta();
175 phi = uArray->GetUnitPhi();
176 theta = EtaToTheta(eta);
177 en = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta));
178 px = TMath::Cos(phi)*pt;
179 py = TMath::Sin(phi)*pt;
180 pz = en*TMath::TanH(eta);
181
182 if (pttrackcut > TMath::Sqrt(px*px+py*py)) continue;
183
184 if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
8838ab7a 185
73faae2f 186 fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object
187 inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm
188 inputParticles.push_back(inputPart); // back of the input_particles vector
189 ipart++;
190 }
191 } // End loop on UnitArray
192 }
8838ab7a 193
194//******************************** CHOICE OF JET AREA
195// Here we determine jets area for subtracting background later
196// For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
197
198 Double_t ghostEtamax = header->GetGhostEtaMax(); // maximum eta in which a ghost can be generated
199 Double_t ghostArea = header->GetGhostArea(); // area of a ghost
200 Int_t activeAreaRepeats = header->GetActiveAreaRepeats(); // do we repeat area calculation?
201 Double_t gridScatter = header->GetGridScatter(); // fractional random fluctuations of the position of the ghosts on the y-phi grid
202 Double_t ktScatter = header->GetKtScatter(); // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid
203 Double_t meanGhostKt = header->GetMeanGhostKt(); // average transverse momentum of the ghosts.
204
205 Double_t areaTypeNumber = header->GetAreaTypeNumber(); // the number determines jet area type
206 fastjet::AreaType areaType = fastjet::active_area;
207 if (areaTypeNumber == 1) areaType = fastjet::active_area;
208 if (areaTypeNumber == 2) areaType = fastjet::active_area_explicit_ghosts;
209 if (areaTypeNumber == 3) areaType = fastjet::one_ghost_passive_area;
210 if (areaTypeNumber == 4) areaType = fastjet::passive_area;
211 if (areaTypeNumber == 5) areaType = fastjet::voronoi_area;
212
213 fastjet::AreaDefinition areaDef;
214
215 if (areaTypeNumber < 5)
2551af9d 216 {
8838ab7a 217 fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
218 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
2551af9d 219 }
220
8838ab7a 221 if (areaTypeNumber == 5)
2551af9d 222 {
8838ab7a 223 Double_t effectiveRFact = header->GetEffectiveRFact();
224 fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
225 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
2551af9d 226 }
227
8838ab7a 228//********************************
229//********************************
2551af9d 230
8838ab7a 231 Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
2551af9d 232
8838ab7a 233//********************************
234//********************************
2551af9d 235
8838ab7a 236//********************************
237//******************************** BG SUBTRACTION
238 if (bgMode == 1)// BG subtraction
73faae2f 239 {
8838ab7a 240 //******************************** JETS FINDING AND EXTRACTION
241 fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
73faae2f 242
8838ab7a 243 // Here we extract inclusive jets with pt > ptmin, sorted by pt
244 Double_t ptMin = header->GetMinJetPt();
245 vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
246 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
247
248 //***************************** BACKGROUND SUBTRACTION
249
250 // Set the rapidity-azimuth range within which to study background
251 Double_t rapMin = header->GetRapMin();
252 Double_t rapMax = header->GetRapMax();
253 Double_t phiMin = header->GetPhiMin();
254 Double_t phiMax = header->GetPhiMax();
255 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
256
73faae2f 257 // As SISCone gives too small areas to estimate background,
258 // we have to choose between kT and Cambridge/Aachen algorithms
259 // to estimate mean background rho
8838ab7a 260 fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
73faae2f 261 Int_t algo = header->GetBGAlgorithm();
8838ab7a 262 if (algo == 0) algorithm = fastjet::kt_algorithm;
263 if (algo == 1) algorithm = fastjet::cambridge_algorithm;
73faae2f 264
265 // Rho estimation :
266 // Radius used to calculate rho, can be different from the one used to reconstruct jets
267 Double_t RRho = header->GetRForRho();
268 fastjet::JetDefinition jetDefinitionForRhoEstimation(algorithm, RRho);
269 fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefinitionForRhoEstimation, areaDef);
270 vector<fastjet::PseudoJet> inclusiveJetsForRhoEstimation = sorted_by_pt(csForRho.inclusive_jets());
271
272 // Number of hard jets not to be used to estimate rho
273 Int_t NHardJets = header->GetNumberOfJetsToErase();
274 if (debug) cout<<"Number of jets before subtraction : "<<inclusiveJetsForRhoEstimation.size()<<endl;
275 if (debug) cout<<"Number of jets not to count into rho estimation : "<<NHardJets<<endl;
276 Bool_t rhoOnJet = 0;
277 if (inclusiveJetsForRhoEstimation.size() == 1 && inclusiveJetsForRhoEstimation[0].perp() > 5) rhoOnJet = 1;//if there's only one jet (not bg kind), rho = 0
6f2a8880 278 if (inclusiveJetsForRhoEstimation.size() <= (UInt_t)NHardJets && inclusiveJetsForRhoEstimation.size() != 1 && inclusiveJetsForRhoEstimation.size() != 0)
73faae2f 279 {
280 // inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin(),inclusiveJetsForRhoEstimation.begin()+1);
281 inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin());
282 cout<<"Only the 1st (hardest) jet of the event hasn't been taken into account for rho estimation"<<endl;
283 }
6f2a8880 284 if (inclusiveJetsForRhoEstimation.size() > (UInt_t)NHardJets) inclusiveJetsForRhoEstimation.erase(inclusiveJetsForRhoEstimation.begin(), inclusiveJetsForRhoEstimation.begin()+NHardJets);
73faae2f 285
286 // Estimation of rho and fluctuations sigma
287
288 //1st method
289 Double_t rho = 0;
290 Double_t sigma = 0;
291 Double_t meanarea = 0;
292 // 3rd argument : 1 = use area 4 vector rather than simple area
293 // last argument : 0 = in case of explicit ghosts use
294 if (inclusiveJetsForRhoEstimation.size() != 0) csForRho.get_median_rho_and_sigma(inclusiveJetsForRhoEstimation, range, 1, rho, sigma, meanarea, 1); // this gives the fluctuations also
295
296 if (rhoOnJet) rho = 0;
297
298 if(debug) cout<<"rho = "<<rho<<endl;
8838ab7a 299
300 // Vector of corrected jets
73faae2f 301 vector<fastjet::PseudoJet> corrJets = sorted_by_pt(clustSeq.subtracted_jets(rho, ptMin));
8838ab7a 302
303 //***************************** JETS DISPLAY
304
305 for (size_t j = 0; j < jets.size(); j++)
306 {
73faae2f 307// // If the jet is only made of ghosts, continue.
308// if (clustSeq.is_pure_ghost(corrJets[j]) == 1) continue;
8838ab7a 309
310 // If the correction is > jet energy px = py = pz = e = 0
311 if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue;
312
73faae2f 313 if(debug)
314 {
315 cout<<"********************************** Reconstructed jet(s) (not corrected)"<<endl;
9b72221b 316 cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<" area : "<<clustSeq.area(jets[j])<<endl;
9b72221b 317 cout<<"e = "<<jets[j].E()<<endl;
318
319 cout<<"********************************** Corrected jet(s)"<<endl;
320 cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<corrJets[j].rap()<<" Phi : "<<corrJets[j].phi()<<" pT : "<<corrJets[j].perp()<<endl;
9b72221b 321 cout<<"e = "<<corrJets[j].E()<<endl;
322 }
8838ab7a 323 // Go to write AOD info
324 AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E());
325 if(debug) aodjet.Print("");
326 AddJet(aodjet);
327 }
328 }
2551af9d 329
8838ab7a 330//********************************
331//******************************** NO BG SUBTRACTION
2551af9d 332
8838ab7a 333 if (bgMode == 0)// No BG subtraction
2551af9d 334 {
8838ab7a 335 //******************************** JETS FINDING AND EXTRACTION
f92ba4be 336 fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
8838ab7a 337 // Here we extract inclusive jets with pt > ptmin, sorted by pt
338 Double_t ptMin = header->GetMinJetPt();
339 vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
340 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
f92ba4be 341
8838ab7a 342 //***************************** JETS DISPLAY
343
344 for (size_t k = 0; k < jets.size(); k++)
345 {
f92ba4be 346 if(debug)
347 {
73faae2f 348 cout<<"********************************** Reconstructed jet(s) (not corrected)"<<endl;
9b72221b 349 cout<<"Jet number "<<k+1<<" : "<<"Rapidity : "<<jets[k].rap()<<" Phi : "<<jets[k].phi()<<" pT : "<<jets[k].perp()<<endl;
9b72221b 350 cout<<"e = "<<jets[k].E()<<endl;
351 }
f92ba4be 352
8838ab7a 353 // Go to write AOD info
f92ba4be 354 Double_t area = clustSeq.area(jets[k]);
8838ab7a 355 AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
f92ba4be 356 aodjet.SetEffArea(area,0);
8838ab7a 357 if(debug) aodjet.Print("");
358 AddJet(aodjet);
359 }
2551af9d 360 }
95345ec0 361
362 delete plugin;
363
2551af9d 364}
365
366//____________________________________________________________________________
367
8838ab7a 368Float_t AliSISConeJetFinder::EtaToTheta(Float_t arg)
369{
370 // return (180./TMath::Pi())*2.*atan(exp(-arg));
371 return 2.*atan(exp(-arg));
372
373
374}
375
376//____________________________________________________________________________
377
378void AliSISConeJetFinder::InitTask(TChain *tree)
379{
380
381 printf("SISCone jet finder initialization ******************");
382 fReader->CreateTasks(tree);
383
384}
385
386
387//____________________________________________________________________________
388
446dbc09 389void AliSISConeJetFinder::WriteJHeaderToFile() const
2551af9d 390{
391 fHeader->Write();
392}
8838ab7a 393
394//____________________________________________________________________________
395