limit printouts
[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
87 AliSISConeJetHeader *header = (AliSISConeJetHeader*)fHeader;
8838ab7a 88 Bool_t debug = header->GetDebug(); // debug option
89 Int_t fOpt = fReader->GetReaderHeader()->GetDetector();
2551af9d 90
91 // Check if we are reading AOD jets
92 TRefArray *refs = 0;
93 Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
94 if (fromAod) { refs = fReader->GetReferences(); }
95
96
97//******************************** SISCONE PLUGIN CONFIGURATION
98// Here we look for SISCone parameters in the header and we define our plugin.
99
100 Double_t coneRadius = header->GetConeRadius(); // cone radius
101 Double_t overlapThreshold = header->GetOverlapThreshold(); // overlap parameter
102 Int_t nPassMax = header->GetNPassMax(); // maximum number of passes
103 Double_t ptProtoJetMin = header->GetPtProtojetMin(); // pT min of protojets
104 Double_t caching = header->GetCaching(); // do we record found cones for this set of data?
105
8838ab7a 106// if (header->GetSplitMergeScale() == 0) fastjet::SISConePlugin::SplitMergeScale splitMergeScale = fastjet::SISConePlugin::SM_pttilde; // There's only one split merge scale
107// Double_t splitMergeStoppingScale = header->GetSplitMergeStoppingScale(); // Additional cut on pt_tilde of protojets
2551af9d 108
109 fastjet::JetDefinition::Plugin * plugin;
110 plugin = new fastjet::SISConePlugin(coneRadius, overlapThreshold, nPassMax, ptProtoJetMin, caching);
111
8838ab7a 112 vector<fastjet::PseudoJet> inputParticles;
2551af9d 113//******************************** READING OF INPUT PARTICLES
114// 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.
115
8838ab7a 116 if(fOpt==0)
117 {
118 TClonesArray *lvArray = fReader->GetMomentumArray();
119 Int_t nIn = lvArray->GetEntries();
120
121 // We check if lvArray is ok
122 if(lvArray == 0)
123 {
124 cout << "Could not get the momentum array" << endl;
125 return;
126 }
127
128 if(nIn == 0)// nIn = Number of particles in the event
129 {
130 if (debug) cout << "entries = 0 ; Event empty !!!" << endl ;
131 return;
132 }
133
8838ab7a 134 Float_t px,py,pz,en;
135
136 // Load input vectors
137 for(Int_t i = 0; i < nIn; i++)
138 {
139 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
140 px = lv->Px();
141 py = lv->Py();
142 pz = lv->Pz();
143 en = lv->Energy();
144
145 fastjet::PseudoJet inputPart(px,py,pz,en);
146 inputPart.set_user_index(i);
147 inputParticles.push_back(inputPart);
148
149 }
150 }
151 else {
152 TClonesArray* fUnit = fReader->GetUnitArray();
153 if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
8838ab7a 154 Int_t nIn = fUnit->GetEntries();
155 if(nIn == 0) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }
42b0ac89 156 // Information extracted from fUnitArray
8838ab7a 157 // load input vectors and calculate total energy in array
158 Float_t pt,eta,phi,theta,px,py,pz,en;
159 Int_t ipart = 0;
160 for(Int_t i=0; i<nIn; i++)
161 {
162 AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
163
164 if(uArray->GetUnitEnergy()>0.){
165
166 // It is not necessary anymore to cut on particle pt
167 pt = uArray->GetUnitEnergy();
168 eta = uArray->GetUnitEta();
169 phi = uArray->GetUnitPhi();
170 theta = EtaToTheta(eta);
171 en = (TMath::Abs(TMath::Sin(theta)) == 0) ? pt : pt/TMath::Abs(TMath::Sin(theta));
172 px = TMath::Cos(phi)*pt;
173 py = TMath::Sin(phi)*pt;
174 pz = en*TMath::TanH(eta);
175 if(debug) cout << "pt: " << pt << ", eta: " << eta << ", phi: " << phi << ", en: " << en << ", px: " << px << ", py: " << py << ", pz: " << pz << endl;
176
177 fastjet::PseudoJet inputPart(px,py,pz,en); // create PseudoJet object
178 inputPart.set_user_index(ipart); //label the particle into Fastjet algortihm
179 inputParticles.push_back(inputPart); // back of the input_particles vector
180 ipart++;
181 }
182 } // End loop on UnitArray
183 }
184
185//******************************** CHOICE OF JET AREA
186// Here we determine jets area for subtracting background later
187// For more informations about jet areas see : The Catchment Area of Jets M. Cacciari, G. Salam and G. Soyez
188
189 Double_t ghostEtamax = header->GetGhostEtaMax(); // maximum eta in which a ghost can be generated
190 Double_t ghostArea = header->GetGhostArea(); // area of a ghost
191 Int_t activeAreaRepeats = header->GetActiveAreaRepeats(); // do we repeat area calculation?
192 Double_t gridScatter = header->GetGridScatter(); // fractional random fluctuations of the position of the ghosts on the y-phi grid
193 Double_t ktScatter = header->GetKtScatter(); // fractional random fluctuations of the tranverse momentum of the ghosts on the y-phi grid
194 Double_t meanGhostKt = header->GetMeanGhostKt(); // average transverse momentum of the ghosts.
195
196 Double_t areaTypeNumber = header->GetAreaTypeNumber(); // the number determines jet area type
197 fastjet::AreaType areaType = fastjet::active_area;
198 if (areaTypeNumber == 1) areaType = fastjet::active_area;
199 if (areaTypeNumber == 2) areaType = fastjet::active_area_explicit_ghosts;
200 if (areaTypeNumber == 3) areaType = fastjet::one_ghost_passive_area;
201 if (areaTypeNumber == 4) areaType = fastjet::passive_area;
202 if (areaTypeNumber == 5) areaType = fastjet::voronoi_area;
203
204 fastjet::AreaDefinition areaDef;
205
206 if (areaTypeNumber < 5)
2551af9d 207 {
8838ab7a 208 fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea, gridScatter, ktScatter, meanGhostKt);
209 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
2551af9d 210 }
211
8838ab7a 212 if (areaTypeNumber == 5)
2551af9d 213 {
8838ab7a 214 Double_t effectiveRFact = header->GetEffectiveRFact();
215 fastjet::VoronoiAreaSpec ghostSpec(effectiveRFact);
216 areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
2551af9d 217 }
218
8838ab7a 219//********************************
220//********************************
2551af9d 221
8838ab7a 222 Bool_t bgMode = header->GetBGMode();// Here one choose to subtract BG or not
2551af9d 223
8838ab7a 224//********************************
225//********************************
2551af9d 226
8838ab7a 227//********************************
228//******************************** BG SUBTRACTION
229 if (bgMode == 1)// BG subtraction
230 {
231 //******************************** JETS FINDING AND EXTRACTION
232 fastjet::ClusterSequenceArea clustSeq(inputParticles, plugin, areaDef);
233 // Here we extract inclusive jets with pt > ptmin, sorted by pt
234 Double_t ptMin = header->GetMinJetPt();
235 vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
236 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
237
238 //***************************** BACKGROUND SUBTRACTION
239
240 // Set the rapidity-azimuth range within which to study background
241 Double_t rapMin = header->GetRapMin();
242 Double_t rapMax = header->GetRapMax();
243 Double_t phiMin = header->GetPhiMin();
244 Double_t phiMax = header->GetPhiMax();
245 fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
246
247 // rho calculus from Cambridge/Aachen clustering (not from SISCONE as it gives too small areas)
248 Int_t algo = header->GetBGAlgorithm();
249 fastjet::JetAlgorithm algorithm = fastjet::kt_algorithm;
250 if (algo == 0) algorithm = fastjet::kt_algorithm;
251 if (algo == 1) algorithm = fastjet::cambridge_algorithm;
252 fastjet::JetDefinition jetDefForRho(algorithm, 0.5);
253 fastjet::ClusterSequenceArea csForRho(inputParticles, jetDefForRho, areaDef);
254 Double_t rho = csForRho.median_pt_per_unit_area_4vector(range);
255 cout<<"rho = "<<rho<<endl;
256
257 // Vector of corrected jets
258 vector<fastjet::PseudoJet> corrJets = clustSeq.subtracted_jets(rho, ptMin);
259
260 //***************************** JETS DISPLAY
261
262 for (size_t j = 0; j < jets.size(); j++)
263 {
264 // If the jet is only made of ghosts, continue.
265 if (clustSeq.is_pure_ghost(jets[j]) == 1) continue;
266
267 // If the correction is > jet energy px = py = pz = e = 0
268 if (corrJets[j].px() == 0 && corrJets[j].py() == 0 && corrJets[j].pz() == 0 && corrJets[j].E() == 0) continue;
269
9b72221b 270 if(debug){
271 cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
272 cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<jets[j].rap()<<" Phi : "<<jets[j].phi()<<" pT : "<<jets[j].perp()<<" area : "<<clustSeq.area(jets[j])<<endl;
273 // cout<<"px = "<<jets[j].px()<<endl;
274 // cout<<"py = "<<jets[j].py()<<endl;
275 // cout<<"pz = "<<jets[j].pz()<<endl;
276 cout<<"e = "<<jets[j].E()<<endl;
277
278 cout<<"********************************** Corrected jet(s)"<<endl;
279 cout<<"Jet number "<<j+1<<" : "<<"Rapidity : "<<corrJets[j].rap()<<" Phi : "<<corrJets[j].phi()<<" pT : "<<corrJets[j].perp()<<endl;
280 // cout<<"px = "<<corrJets[j].px()<<endl;
281 // cout<<"py = "<<corrJets[j].py()<<endl;
282 // cout<<"pz = "<<corrJets[j].pz()<<endl;
283 cout<<"e = "<<corrJets[j].E()<<endl;
284 }
8838ab7a 285 // Go to write AOD info
286 AliAODJet aodjet (corrJets[j].px(), corrJets[j].py(), corrJets[j].pz(), corrJets[j].E());
287 if(debug) aodjet.Print("");
288 AddJet(aodjet);
289 }
290 }
2551af9d 291
8838ab7a 292//********************************
293//******************************** NO BG SUBTRACTION
2551af9d 294
8838ab7a 295 if (bgMode == 0)// No BG subtraction
2551af9d 296 {
8838ab7a 297 //******************************** JETS FINDING AND EXTRACTION
298 fastjet::ClusterSequence clustSeq(inputParticles, plugin);
299 // Here we extract inclusive jets with pt > ptmin, sorted by pt
300 Double_t ptMin = header->GetMinJetPt();
301 vector<fastjet::PseudoJet> inclusiveJets = clustSeq.inclusive_jets(ptMin);
302 vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets);
303
304 //***************************** JETS DISPLAY
305
306 for (size_t k = 0; k < jets.size(); k++)
307 {
9b72221b 308 if(debug){
309 cout<<"********************************** Reconstructed jet(s) (non corrected)"<<endl;
310 cout<<"Jet number "<<k+1<<" : "<<"Rapidity : "<<jets[k].rap()<<" Phi : "<<jets[k].phi()<<" pT : "<<jets[k].perp()<<endl;
311 // cout<<"px = "<<jets[k].px()<<endl;
312 // cout<<"py = "<<jets[k].py()<<endl;
313 // cout<<"pz = "<<jets[k].pz()<<endl;
314 cout<<"e = "<<jets[k].E()<<endl;
315 }
8838ab7a 316 // Go to write AOD info
317 AliAODJet aodjet (jets[k].px(), jets[k].py(), jets[k].pz(), jets[k].E());
318 if(debug) aodjet.Print("");
319 AddJet(aodjet);
320 }
2551af9d 321 }
322}
323
8838ab7a 324//____________________________________________________________________________
325
326Float_t AliSISConeJetFinder::EtaToTheta(Float_t arg)
327{
328 // return (180./TMath::Pi())*2.*atan(exp(-arg));
329 return 2.*atan(exp(-arg));
330
331
332}
333
334//____________________________________________________________________________
335
336void AliSISConeJetFinder::InitTask(TChain *tree)
337{
338
339 printf("SISCone jet finder initialization ******************");
340 fReader->CreateTasks(tree);
341
342}
343
344
2551af9d 345//____________________________________________________________________________
346
347void AliSISConeJetFinder::WriteJHeaderToFile()
348{
349 fHeader->Write();
350}
8838ab7a 351
352//____________________________________________________________________________
353