]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliFastJetFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
In AliEveEventManager::AssertGeometry() remove a hack to bypass
[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
16//---------------------------------------------------------------------
17// FastJet finder
18// interface to FastJet algorithm
19// Author: Rafael.Diaz.Valdes@cern.ch
20// kt using NlnN
21//---------------------------------------------------------------------
22
23#include <Riostream.h>
24#include <TLorentzVector.h>
25#include <TFile.h>
26#include <TH1F.h>
27#include <TH2F.h>
28#include <TArrayF.h>
29#include <TRandom.h>
30#include <TClonesArray.h>
31
32#include "AliFastJetFinder.h"
33#include "AliFastJetHeader.h"
34#include "AliJetReaderHeader.h"
35#include "AliJetReader.h"
36#include "AliJet.h"
37//for FastJet finder
38#include "FjPseudoJet.hh"
39#include "FjClusterSequence.hh"
40
41
42ClassImp(AliFastJetFinder);
43
44
45////////////////////////////////////////////////////////////////////////
46
47AliFastJetFinder::AliFastJetFinder():
48 AliJetFinder(),
49 fHeader(0x0),
50 fLego(0x0),
51 fLegoSignal(0x0)
52{
53 // Constructor
54}
55
56////////////////////////////////////////////////////////////////////////
57
58AliFastJetFinder::~AliFastJetFinder()
59
60{
61 // destructor
62}
63
64////////////////////////////////////////////////////////////////////////
65
66
67void AliFastJetFinder::FindJets()
68
69{
70 //create cells and jets array
71 // 1) transform input to pt,eta,phi plus lego
72 TClonesArray *lvArray = fReader->GetMomentumArray();
73 Int_t nIn = lvArray->GetEntries();
74 if (nIn == 0) return;
75
76 // local arrays for particles input
77 Float_t* enT = new Float_t[nIn];
78 Float_t* ptT = new Float_t[nIn];
79 Float_t* etaT = new Float_t[nIn];
80 Float_t* phiT = new Float_t[nIn];
81 Int_t* injet = new Int_t[nIn];
82
83
84 //total energy in array
85 Float_t EtTotal = 0.0;
86 Float_t meanptCell = 0.0;
87 Float_t sqptCell = 0.0;
88
89
90 // load input vectors in fLego
91 for (Int_t i = 0; i < nIn; i++){
92 TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
93 enT[i] = lv->Energy();
94 ptT[i] = lv->Pt();
95 etaT[i] = lv->Eta();
96 phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
97 if (fReader->GetCutFlag(i) == 1){
98 fLego->Fill(etaT[i], phiT[i], ptT[i]);
99 if(fReader->GetSignalFlag(i) == 1)
100 fLegoSignal->Fill(etaT[i], phiT[i], ptT[i]);
101 EtTotal= EtTotal+ptT[i];
102 }
103 }
104 fJets->SetNinput(nIn);
105
106 // add soft background fixed
107 Int_t nsoft = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
108 Float_t* ptRndm = new Float_t[nsoft];
109 if(fHeader->AddSoftBackg()){
110 gRandom->RndmArray(nsoft,ptRndm);
111 for(Int_t isoft = 0; isoft < nsoft; isoft++){
112 Float_t ptsoft = 0.005*ptRndm[isoft];
113 EtTotal= EtTotal+ptsoft;
114 }
115 }
116
117 if(EtTotal == 0){
118 delete enT;
119 delete ptT;
120 delete etaT;
121 delete phiT;
122 return;
123 }
124
125 //cell array
126 Float_t* etCell = new Float_t[90000]; //! Cell Energy // check enough space! *to be done*
127 Float_t* etaCell = new Float_t[90000]; //! Cell eta
128 Float_t* phiCell = new Float_t[90000]; //! Cell phi
129 Int_t* jetflagCell = new Int_t[90000]; //! Cell flag for jets
130 Float_t* etsigCell = new Float_t[90000]; // signal in this cell
131
132 //jets array
133 Float_t* etaJet = new Float_t[200];
134 Float_t* phiJet = new Float_t[200];
135 Float_t* etJet = new Float_t[200];
136 Float_t* etsigJet = new Float_t[200]; //signal energy
137 Float_t* etallJet = new Float_t[200]; //total energy in jet area
138 Int_t* ncellsJet = new Int_t[200];
139 memset(etaJet,0,sizeof(Float_t)*200);
140 memset(phiJet,0,sizeof(Float_t)*200);
141 memset(etJet,0,sizeof(Float_t)*200);
142 memset(etsigJet,0,sizeof(Float_t)*200);
143 memset(etallJet,0,sizeof(Float_t)*200);
144 memset(ncellsJet,0,sizeof(Int_t)*200);
145
146
147
148 // load cells arrays
149 Int_t nCell = 0;
150 TAxis* xaxis = fLego->GetXaxis();
151 TAxis* yaxis = fLego->GetYaxis();
152 Float_t e = 0.0; Float_t esig = 0.0;
153
154 for (Int_t ib = 1; ib <= fHeader->GetLegoNbinEta(); ib++) {
155 for (Int_t jb = 1; jb <= fHeader->GetLegoNbinPhi(); jb++) {
156 e = fLego->GetBinContent(ib,jb);
157 if (e < 0.0) continue; // don't include this cells
158 Float_t eta = xaxis->GetBinCenter(ib);
159 Float_t phi = yaxis->GetBinCenter(jb);
160 if(fHeader->AddSoftBackg())
161 etCell[nCell] = e + 0.005*ptRndm[nCell];
162 else
163 etCell[nCell] = e;
164 sqptCell = sqptCell + etCell[nCell]*etCell[nCell]; // xi^2 ////////
165 etaCell[nCell] = eta;
166 phiCell[nCell] = phi;
167 jetflagCell[nCell] = -1; //default
168 esig = fLegoSignal->GetBinContent(ib,jb);
169 if(esig > 0.0)
170 etsigCell[nCell] = esig;
171 else
172 etsigCell[nCell] = 0.0;
173 nCell++;
174 }
175 }
176
177 meanptCell = EtTotal/(Float_t)nCell;
178 sqptCell = sqptCell/(Float_t)nCell;
179
180 Int_t nJets = 0;
181 //call to FastJet Algorithm
182 RunAlgorithm(nJets,etJet,etaJet,phiJet,etsigJet,etallJet,ncellsJet,
183 nCell,etCell,etaCell,phiCell,etsigCell,jetflagCell);
184
185
186 //subtract background
187 SubtractBackg(nCell,jetflagCell,etCell,
188 nJets,etJet,etallJet,ncellsJet,
189 meanptCell,sqptCell,EtTotal);
190
191
192 // add jets to list
193 Int_t* index = new Int_t[nJets];
194 Int_t nj = 0;
195 for(Int_t kj=0; kj<nJets; kj++){
196 if ((etaJet[kj] > (fHeader->GetJetEtaMax())) ||
197 (etaJet[kj] < (fHeader->GetJetEtaMin())) ||
198 (etJet[kj] < fHeader->GetMinJetEt())) continue; // acceptance eta range and etmin
199 Float_t px, py,pz,en; // convert to 4-vector
200 px = etJet[kj] * TMath::Cos(phiJet[kj]);
201 py = etJet[kj] * TMath::Sin(phiJet[kj]);
202 pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj])));
203 en = TMath::Sqrt(px * px + py * py + pz * pz);
204 fJets->AddJet(px, py, pz, en);
205 index[nj] = kj;
206 nj++;
207 }
208
209 //add signal percentage and total signal in AliJets for analysis tool
210 Float_t* percentage = new Float_t[nj];
211 Int_t* ncells = new Int_t[nj];
212 Int_t* mult = new Int_t[nj];
213
214 for(Int_t i = 0; i< nj; i++){
215 percentage[i] = etsigJet[index[i]]/etJet[index[i]];
216 ncells[i] = ncellsJet[index[i]];
217 }
218
219 //reorder injet flags
220 for(Int_t ipar = 0; ipar < nIn; ipar++){
221 Float_t injetflag =0;
222 Int_t iparCell = fLego->FindBin(etaT[ipar], phiT[ipar]);
223 injet[ipar] = jetflagCell[iparCell];
224 for(Int_t js = 0; js < nj; js++){
225 if(injet[ipar] == index[js]){
226 injet[ipar] = js; // set the new jet id value
227 mult[js]++; // add multiplicity in jet js
228 injetflag = 1;
229 break;
230 }
231 }
232 if(injetflag == 0) injet[ipar] = -1; // set default value
233 }
234
235
236 fJets->SetNCells(ncells);
237 fJets->SetPtFromSignal(percentage);
238 fJets->SetMultiplicities(mult);
239 fJets->SetInJet(injet);
240 fJets->SetEtaIn(etaT);
241 fJets->SetPhiIn(phiT);
242 fJets->SetPtIn(ptT);
243 fJets->SetEtAvg(meanptCell);
244
245 //delete
246 delete enT;
247 delete ptT;
248 delete etaT;
249 delete phiT;
250 delete injet;
251 //cells
252 delete etCell;
253 delete etaCell;
254 delete phiCell;
255 delete jetflagCell;
256 delete etsigCell;
257 //jets
258 delete etaJet;
259 delete phiJet;
260 delete etJet;
261 delete etsigJet;
262 delete etallJet;
263 delete ncellsJet;
264
265 delete index;
266 delete percentage;
267 delete ncells;
268 delete mult;
269 delete ptRndm;
270
271}
272
273////////////////////////////////////////////////////////////////////////
274void AliFastJetFinder::RunAlgorithm(Int_t& nJets,Float_t* etJet,Float_t* etaJet,Float_t* phiJet,
275 Float_t* etsigJet, Float_t* etallJet, Int_t* ncellsJet,
276 Int_t& nCell,Float_t* etCell,Float_t* etaCell,Float_t* phiCell,
277 Float_t* etsigCell, Int_t* jetflagCell)
278{
279 //FastJet objects
280 vector<FjPseudoJet> input_cells; // create a vector
281 for (Int_t i = 0; i < nCell; i++){
282 if(etCell[i] == 0.0) continue; // not include cell empty
283 Double_t px, py,pz,en; // convert to 4-vector
284 px = etCell[i]*TMath::Cos(phiCell[i]);
285 py = etCell[i]*TMath::Sin(phiCell[i]);
286 pz = etCell[i]/TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaCell[i])));
287 en = TMath::Sqrt(px * px + py * py + pz * pz);
288 FjPseudoJet input_cell(px,py,pz,en); // create FjPseudoJet object
289 input_cell.set_user_index(i); //label the cell into Fastjet algortihm
290 //push FjPseudoJet of (px,py,pz,en) onto back of the input_cells
291 input_cells.push_back(input_cell);
292 }
293
294 //run the jet clustering with option R=1.0 and strategy= Best
295 Double_t Rparam = fHeader->GetRadius(); // default 1.0;
296 FjClusterSequence clust_seq(input_cells,Rparam);
297
298
299 //vector to get clusters
300 vector<FjPseudoJet> clusters;
301
302 ///////////////////////////////////////////////////////////////////////////
303 //extract the inclusive jets with pt> ptmin sorted by pt
304 //Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
305 // Double_t ptbgRc = EtBackgT*(Rparam*Rparam*TMath::Pi()/areaT);
306 // Double_t ptbgRcfluct = dEtTotal*Rparam*TMath::Sqrt(TMath::Pi()/areaT);
307 // Double_t ptmin = ptbgRc + ptbgRcfluct;
308 clusters = clust_seq.inclusive_jets(0);
309 //////////////////////////////////////////////////////////////////////////
310
311 /////////////////////////////////////////////////////////////////////////
312 //extract the exclusive jets with dcut = 25 GeV**2 and sort them in order
313 //of increasing pt
314 //Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
315 //Double_t ptbgRc = EtBackgT*(Rparam*Rparam*TMath::Pi()/areaT);
316 //Double_t ptbgRcfluct = dEtTotal*Rparam*TMath::Sqrt(TMath::Pi()/areaT);
317 //Double_t ptmin = ptbgRc + ptbgRcfluct;
318 //Double_t ktbackg = (fHeader->GetKfactor())*ptmin;
319 //Double_t dcut = ktbackg*ktbackg;
320 //clusters = sorted_by_pt(clust_seq.exclusive_jets(dcut));
321 //clusters = sorted_by_pt(clust_seq.exclusive_jets(5));
322 /////////////////////////////////////////////////////////////////////////
323
324
325 //cout << " ktClusters " << clusters.size() << endl;
326 nJets = (Int_t)clusters.size();
327 ////////////////////////////////////////////////////////////////////////
328 // add all clusters to jet arrays
329 for(UInt_t ij = 0; ij < clusters.size(); ij++){
330 //constituents
331 vector<FjPseudoJet> constituents = clust_seq.constituents(clusters[ij]);
332 //fill jet array info
333 ncellsJet[ij] = (Int_t)constituents.size();
334 phiJet[ij] = clusters[ij].phi();
335 Float_t angle = TMath::ATan(clusters[ij].perp()/clusters[ij].pz());
336 angle = ((angle < 0) ? angle + TMath::Pi() : angle);
337 etaJet[ij] = - TMath::Log(TMath::Tan(angle/2.0));
338 etJet[ij] = clusters[ij].perp();
339 //get constituents cells
340 for(UInt_t jc = 0; jc < constituents.size(); jc++){ // loop for all cells in ij cluster
341 Int_t jcell = constituents[jc].user_index();
342 jetflagCell[jcell] = ij; //flag this cell for jet
343 etsigJet[ij] = etsigJet[ij] + etsigCell[jcell]; // add signal of this cell
344 etallJet[ij] = etallJet[ij] + etCell[jcell]; // add total of this cell
345 }
346 }
347
348}
349////////////////////////////////////////////////////////////////////////
350void AliFastJetFinder::SubtractBackg(Int_t& nCell, Int_t* jetflagCell, Float_t* etCell,
351 Int_t& nJets, Float_t* etJet, Float_t* etallJet, Int_t* ncellsJet,
352 Float_t& meanptCell, Float_t& sqptCell, Float_t& etBackg)
353{
354 // simplest method: subtract background from external region to jets
355
356 //tmp array to flag jets
357 Int_t flagJet[200];
358 //tmp array to flag jet-cell status
359 Int_t tmpjetflagCell[90000];
360
361 Float_t etBackgOld = 0;
362 Float_t prec = fHeader->GetPrecBg();
363 Float_t bgprec = 1;
364
365 while(bgprec > prec){
366 //clear tmpjetflagCell
367 memset(tmpjetflagCell,-1,sizeof(Int_t)*90000); // init with -1 (all cells are background)
368 //clear flagjet
369 memset(flagJet,0,sizeof(Int_t)*200); // init with 0 (no flag jets)
370 // select clusters > meantmpCell
371 for(Int_t i = 0; i < nJets; i++){
372 Float_t iptcell = etallJet[i]/(Float_t)ncellsJet[i];
373 if(iptcell < meanptCell) continue; // cluster not selected
374 // convert tmp cell background to jet cell
375 for(Int_t ic = 0; ic < nCell; ic++){ //loop for all cells
376 if(jetflagCell[ic] != i) continue; // other cells
377 tmpjetflagCell[ic] = i; // convert to a jet cell
378 }
379 //load total energy in cluster
380 etJet[i] = etallJet[i];
381 flagJet[i] = 1; // flag jet
382 }
383 //subtract background
384 for(Int_t j = 0; j < nCell; j++){ // loop for all cells
385 Int_t idxjet = tmpjetflagCell[j];
386 if(idxjet == -1) continue; // background cell
387 if(idxjet == -2) continue; // background protected cell
388 etJet[idxjet] = etJet[idxjet] - meanptCell;
389 }
390 // evaluate background fluctuations (rms value)
391 Float_t rmsptCell = TMath::Sqrt(sqptCell - meanptCell*meanptCell);
392 //fake jets
393 for(Int_t k = 0; k < nJets; k++){
394 if(flagJet[k] != 1) continue; // only flaged jets
395 //if(etJet[k] > fHeader->GetMinJetEt()) continue; // jet ok!!
396 if(etJet[k] > rmsptCell*ncellsJet[k]) continue; // jet ok!!
397 //clear tmpjetflag in cells
398 for(Int_t kc = 0; kc < nCell; kc++){ //loop for all cells
399 if(tmpjetflagCell[kc] != k) continue; // other cells
400 tmpjetflagCell[kc] = -1; // convert to background tmp cell
401 }
402 // clear all previous jet flags
403 etJet[k] = 0;
404 flagJet[k] = 0;
405 }
406 // recalculate background
407 sqptCell = 0;
408 etBackgOld = etBackg;
409 etBackg = 0;
410 Int_t nCellBackg = 0;
411 for(Int_t l = 0; l < nCell; l++){ // loop for all cells
412 if(tmpjetflagCell[l] != -1) continue; // cell included in some jet or protected
413 nCellBackg++;
414 etBackg = etBackg + etCell[l]; //add cell to background
415 //calc sqptCell
416 sqptCell = sqptCell + etCell[l]*etCell[l];
417 }
418 if(nCellBackg){
419 meanptCell = etBackg/(Float_t)nCellBackg; // new pt cell mean value
420 sqptCell = sqptCell/(Float_t)nCellBackg;
421 }else{
422 meanptCell = 0;
423 sqptCell = 0;
424 }
425 // evaluate presicion values
426 if(etBackg)
427 bgprec = (etBackgOld - etBackg)/etBackg;
428 else
429 bgprec = 0;
430 }
431
432 // set etJet 0 for all clusters not flaged in order to
433 for(Int_t m = 0; m < nJets; m++){
434 if(flagJet[m] == 1) continue; // flaged jets
435 etJet[m] = 0; //others clusters
436 }
437
438
439}
440
441////////////////////////////////////////////////////////////////////////
442void AliFastJetFinder::SubtractBackgArea(Int_t& nCell, Int_t* jetflagCell,
443 Float_t* etCell,Int_t&nJets, Float_t* etJet, Float_t* etallJet)
444{
445 // area method: subtract background from external region to jets
446 // using fixed area pi*Rc2
447
448 // n cells contained in a cone Rc
449 Double_t Rc = fHeader->GetRadius();
450 Float_t nCellRc = Rc*Rc*TMath::Pi()/(0.015*0.015); // change in future !!!!
451 //tmp array to flag fake jets
452 Int_t fakeflagJet[100];
453 memset(fakeflagJet,0,sizeof(Int_t)*100); // init with 0 (no fake jets)
454 Int_t njfake = nJets;
455 while(njfake > 0){
456 //calc background per cell
457 Int_t nCellBackg = 0;
458 Float_t EtBackg = 0.0;
459 for(Int_t i = 0; i < nCell; i++){ // loop for all cells
460 if(jetflagCell[i] != -1) continue; // cell included in some jet
461 nCellBackg++;
462 EtBackg = EtBackg + etCell[i]; //add cell to background
463 }
464 //subtract background energy per jet
465 for(Int_t l = 0; l < nJets; l++){
466 if(fakeflagJet[l] == 1) continue; // fake jet
467 etJet[l] = etallJet[l] - nCellRc*EtBackg/(Float_t)nCellBackg;
468 }
469 //fake jets analysis
470 njfake = 0;
471 for(Int_t k = 0; k < nJets; k++){
472 if(fakeflagJet[k] == 1) continue; // fake jet
473 if(etJet[k] < fHeader->GetMinJetEt()){ // a new fake jet
474 //clear jet flag in cells
475 for(Int_t kc = 0; kc < nCell; kc++){ //loop for all cells
476 Int_t kidx = jetflagCell[kc];
477 if(kidx != k) continue; // other cells
478 jetflagCell[kc] = -1; // convert to background cell
479 }
480 fakeflagJet[k] = 1; // mark as a fake jet
481 njfake++; //count fakes in this loop
482 }
483 }
484 }
485}
486
487////////////////////////////////////////////////////////////////////////
488
489
490
491void AliFastJetFinder::Reset()
492{
493 fLego->Reset();
494 fLegoSignal->Reset();
495 fJets->ClearJets();
496
497}
498////////////////////////////////////////////////////////////////////////
499
500void AliFastJetFinder::WriteJHeaderToFile()
501{
502 fOut->cd();
503 fHeader->Write();
504}
505
506////////////////////////////////////////////////////////////////////////
507
508void AliFastJetFinder::Init()
509{
510 // initializes some variables
511 // book lego
512 fLego = new
513 TH2F("legoH","eta-phi",
514 fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
515 fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
516 fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
517 fLegoSignal = new
518 TH2F("legoSignalH","eta-phi signal",
519 fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
520 fHeader->GetLegoEtaMax(), fHeader->GetLegoNbinPhi(),
521 fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
522
523}