]>
Commit | Line | Data |
---|---|---|
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 | ||
d2eb0695 | 53 | #if defined(ENABLE_PLUGIN_SISCONE) || defined(FASTJET_ENABLE_PLUGIN_SISCONE) |
2551af9d | 54 | #include "fastjet/SISConePlugin.hh" |
55 | #endif | |
56 | ||
8838ab7a | 57 | #include<sstream> // needed for internal io |
2551af9d | 58 | #include<vector> |
59 | #include <cmath> | |
60 | ||
61 | using namespace std; | |
62 | ||
8838ab7a | 63 | |
2551af9d | 64 | ClassImp(AliSISConeJetFinder) |
65 | ||
8838ab7a | 66 | |
2551af9d | 67 | //____________________________________________________________________________ |
68 | ||
69 | AliSISConeJetFinder::AliSISConeJetFinder(): | |
70 | AliJetFinder() | |
71 | { | |
72 | // Constructor | |
73 | } | |
74 | ||
75 | //____________________________________________________________________________ | |
76 | ||
77 | AliSISConeJetFinder::~AliSISConeJetFinder() | |
78 | { | |
79 | // destructor | |
80 | } | |
81 | ||
82 | //______________________________________________________________________________ | |
83 | void 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 | ||
8838ab7a | 366 | //____________________________________________________________________________ |
367 | ||
368 | Float_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 | ||
378 | void AliSISConeJetFinder::InitTask(TChain *tree) | |
379 | { | |
380 | ||
381 | printf("SISCone jet finder initialization ******************"); | |
382 | fReader->CreateTasks(tree); | |
383 | ||
384 | } | |
385 | ||
386 | ||
2551af9d | 387 | //____________________________________________________________________________ |
388 | ||
446dbc09 | 389 | void AliSISConeJetFinder::WriteJHeaderToFile() const |
2551af9d | 390 | { |
391 | fHeader->Write(); | |
392 | } | |
8838ab7a | 393 | |
394 | //____________________________________________________________________________ | |
395 |