]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/AliUA1JetFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add the possibility to process the event stat file
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinder.cxx
CommitLineData
d89b8229 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/* $Id */
17
18//---------------------------------------------------------------------
19// UA1 Cone Algorithm Jet finder for jet studies
20// manages the search for jets using charged particle momentum and
21// neutral cell energy information
22// Authors: Rafael.Diaz.Valdes@cern.ch
23// magali.estienne@subatech.in2p3.fr
24// alexandre.shabetai@cern.ch
25// ** 2011
26// Modified accordingly to reader/finder splitting and new handling of neutral information
27// Versions V1 and V2 merged
28//---------------------------------------------------------------------
29
30#include <TH2F.h>
31#include <TMath.h>
32
33#include "AliUA1JetFinder.h"
34#include "AliUA1JetHeaderV1.h"
35#include "AliJetCalTrk.h"
36#include "AliJetBkg.h"
37#include "AliAODJetEventBackground.h"
38#include "AliAODJet.h"
39
40ClassImp(AliUA1JetFinder)
41
42////////////////////////////////////////////////////////////////////////
43
44AliUA1JetFinder::AliUA1JetFinder():
45 AliJetFinder(),
46 fLego(0),
47 fJetBkg(new AliJetBkg())
48{
49 // Default constructor
50}
51
52//-----------------------------------------------------------------------
53AliUA1JetFinder::~AliUA1JetFinder()
54{
55 // Destructor
56 delete fLego;
57 delete fJetBkg;
58
59}
60
61//-----------------------------------------------------------------------
62void AliUA1JetFinder::FindJets()
63{
64 // Used to find jets using charged particle momentum information
65 // & neutral energy from calo cells
66 //
67 // 1) Fill cell map array
68 // 2) calculate total energy and fluctuation level
69 // 3) Run algorithm
70 // 3.1) look centroides in cell map
71 // 3.2) calculate total energy in cones
72 // 3.3) flag as a possible jet
73 // 3.4) reorder cones by energy
74 // 4) subtract backg in accepted jets
75 // 5) fill AliJet list
76
77 // transform input to pt,eta,phi plus lego
78
79 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
80 Int_t nIn = fCalTrkEvent->GetNCalTrkTracks();
81 fDebug = fHeader->GetDebug();
82
83 if (nIn <= 0) return;
84 fJetBkg->SetHeader(fHeader);
85 fJetBkg->SetCalTrkEvent(GetCalTrkEvent());
86 fJetBkg->SetDebug(fDebug);
87 // local arrays for input
88 // ToDo: check memory fragmentation, maybe better to
89 // define them globally and resize as needed
90 // Fragmentation should be worse for low mult...
91 Float_t* ptT = new Float_t[nIn];
92 Float_t* etaT = new Float_t[nIn];
93 Float_t* phiT = new Float_t[nIn];
94 Int_t* injet = new Int_t[nIn];
95 Int_t* injetOk = new Int_t[nIn];
96
97 memset(ptT,0,sizeof(Float_t)*nIn);
98 memset(etaT,0,sizeof(Float_t)*nIn);
99 memset(phiT,0,sizeof(Float_t)*nIn);
100 memset(injet,0,sizeof(Int_t)*nIn);
101 memset(injetOk,-1,sizeof(Int_t)*nIn);
102
103 // load input vectors and calculate total energy in array
104
105 // total energy in array
106 Float_t etbgTotal = 0.;
107 Float_t npart = 0.;
108 Float_t etbg2 = 0.;
109
110 for (Int_t i = 0; i < fCalTrkEvent->GetNCalTrkTracks(); i++){
111 ptT[i] = fCalTrkEvent->GetCalTrkTrack(i)->GetPt();
112 etaT[i] = fCalTrkEvent->GetCalTrkTrack(i)->GetEta();
113 phiT[i] = ((fCalTrkEvent->GetCalTrkTrack(i)->GetPhi() < 0) ? (fCalTrkEvent->GetCalTrkTrack(i)->GetPhi()) + 2 * TMath::Pi() : fCalTrkEvent->GetCalTrkTrack(i)->GetPhi());
114 //fCalTrkEvent->GetCalTrkTrack(i)->Print(Form("%d",i));
115 if (fCalTrkEvent->GetCalTrkTrack(i)->GetCutFlag() != 1) continue;
116 fLego->Fill(etaT[i], phiT[i], ptT[i]);
117 npart += 1;
118 etbgTotal+= ptT[i];
119 etbg2 += ptT[i]*ptT[i];
120 }
121
122 // calculate total energy and fluctuation in map
123 Double_t meanpt = 0.;
124 Double_t ptRMS = 0.;
125 if(npart>0){
126 meanpt = etbgTotal/npart;
127 etbg2 = etbg2/npart;
128 if(etbg2>(meanpt*meanpt)){// prenent NAN, should only happen due to numerical instabilities
129 ptRMS = TMath::Sqrt(etbg2-meanpt*meanpt);
130 }
131 }
132 Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
133
134 // arrays to hold jets
135 Float_t etaJet[kMaxJets];
136 Float_t phiJet[kMaxJets];
137 Float_t etJet[kMaxJets];
138 Float_t etsigJet[kMaxJets]; //signal et in jet
139 Float_t etallJet[kMaxJets]; // total et in jet (tmp variable)
140 Int_t ncellsJet[kMaxJets];
141 Int_t multJetT[kMaxJets];
142 Int_t multJetC[kMaxJets];
143 Int_t multJet[kMaxJets];
144 Float_t *areaJet = new Float_t[kMaxJets];
145 // Used for jet reordering at the end of the jet finding procedure
146 Float_t etaJetOk[kMaxJets];
147 Float_t phiJetOk[kMaxJets];
148 Float_t etJetOk[kMaxJets];
149 Float_t etsigJetOk[kMaxJets]; //signal et in jet
150 Float_t etallJetOk[kMaxJets]; // total et in jet (tmp variable)
151 Int_t ncellsJetOk[kMaxJets];
152 Int_t multJetOk[kMaxJets];
153 Float_t *areaJetOk = new Float_t[kMaxJets];
154 Int_t nJets; // to hold number of jets found by algorithm
155 Int_t nj; // number of jets accepted
156 Float_t prec = header->GetPrecBg();
157 Float_t bgprec = 1;
158
159 while(bgprec > prec){
160 //reset jet arrays in memory
161 memset(etaJet,0,sizeof(Float_t)*kMaxJets);
162 memset(phiJet,0,sizeof(Float_t)*kMaxJets);
163 memset(etJet,0,sizeof(Float_t)*kMaxJets);
164 memset(etallJet,0,sizeof(Float_t)*kMaxJets);
165 memset(etsigJet,0,sizeof(Float_t)*kMaxJets);
166 memset(ncellsJet,0,sizeof(Int_t)*kMaxJets);
167 memset(multJetT,0,sizeof(Int_t)*kMaxJets);
168 memset(multJetC,0,sizeof(Int_t)*kMaxJets);
169 memset(multJet,0,sizeof(Int_t)*kMaxJets);
170 memset(areaJet,0,sizeof(Float_t)*kMaxJets);
171 memset(etaJetOk,0,sizeof(Float_t)*kMaxJets);
172 memset(phiJetOk,0,sizeof(Float_t)*kMaxJets);
173 memset(etJetOk,0,sizeof(Float_t)*kMaxJets);
174 memset(etallJetOk,0,sizeof(Float_t)*kMaxJets);
175 memset(etsigJetOk,0,sizeof(Float_t)*kMaxJets);
176 memset(ncellsJetOk,0,sizeof(Int_t)*kMaxJets);
177 memset(multJetOk,0,sizeof(Int_t)*kMaxJets);
178 memset(areaJetOk,0,sizeof(Float_t)*kMaxJets);
179 nJets = 0;
180 nj = 0;
181 // reset particles-jet array in memory
182 memset(injet,-1,sizeof(Int_t)*nIn);
183 //run cone algorithm finder
184 RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
185 //run background subtraction
186 if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
187 nj = header->GetNAcceptJets();
188 else
189 nj = nJets;
190
191 //subtract background
192 Float_t etbgTotalN = 0.0; //new background
193 Float_t sigmaN = 0.0; //new background
194 if(header->GetBackgMode() == 1) {// standard
195 fJetBkg->SubtractBackg(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);}
196 if(header->GetBackgMode() == 2) //cone
197 fJetBkg->SubtractBackgCone(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
198 if(header->GetBackgMode() == 3) //ratio
199 fJetBkg->SubtractBackgRatio(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
200 if(header->GetBackgMode() == 4) //statistic
201 fJetBkg->SubtractBackgStat(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
202
203 //calc precision
204 if(etbgTotalN != 0.0)
205 bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
206 else
207 bgprec = 0;
208 etbgTotal = etbgTotalN; // update with new background estimation
209
210 } //end while
211
212 // add tracks to the jet if it wasn't yet done
213 if (header->GetBackgMode() == 0){
214 Float_t rc= header->GetRadius();
215 for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
216 for(Int_t ijet=0; ijet<nJets; ijet++){
217 Float_t deta = etaT[jpart] - etaJet[ijet];
218 Float_t dphi = phiT[jpart] - phiJet[ijet];
219 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
220 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
221 Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
222 if(dr <= rc){ // particles inside this cone
223 injet[jpart] = ijet;
224 break;
225 }
226 }// end jets loop
227 } //end particle loop
228 }
229
230 // add jets to list
231 if (fDebug>1) printf("Found %d jets \n", nj);
232
233 // Reorder jets by et in cone
234 // Sort jets by energy
235 Int_t idx[kMaxJets];
236 TMath::Sort(nJets, etJet, idx);
237 for(Int_t p = 0; p < nJets; p++)
238 {
239 etaJetOk[p] = etaJet[idx[p]];
240 phiJetOk[p] = phiJet[idx[p]];
241 etJetOk[p] = etJet[idx[p]];
242 etallJetOk[p] = etJet[idx[p]];
243 etsigJetOk[p] = etsigJet[idx[p]];
244 ncellsJetOk[p] = ncellsJet[idx[p]];
245 multJetOk[p] = multJet[idx[p]];
246 areaJetOk[p] = areaJet[idx[p]];
247 }
248
249 //////////////////////////
250
251 Int_t nTracks = fCalTrkEvent->GetNCalTrkTracks();
252
253 for(Int_t kj=0; kj<nj; kj++)
254 {
255 if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
256 (etaJetOk[kj] < (header->GetJetEtaMin())) ||
257 (phiJetOk[kj] > (header->GetJetPhiMax())) ||
258 (phiJetOk[kj] < (header->GetJetPhiMin())) ||
259 (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
260 Float_t px=-999, py=-999 ,pz=-999 ,en=-999; // convert to 4-vector
261 px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
262 py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
263 pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
264 en = TMath::Sqrt(px * px + py * py + pz * pz);
265 AliAODJet jet(px, py, pz, en);
266
267 // Calc jet area if it wasn't yet done
268 if (header->GetBackgMode() == 0){
269 // calculate the area of the jet
270 Float_t rc= header->GetRadius();
271 areaJetOk[kj] = fJetBkg->CalcJetAreaEtaCut(rc,etaJetOk[kj]);
272 }
273
274 jet.SetEffArea(areaJetOk[kj],0.,0.,0.);
275 jet.SetBgEnergy(etbgTotal,0.);
276 if (fDebug>1) jet.Print(Form("%d",kj));
277
278 for(Int_t jpart = 0; jpart < nTracks; jpart++) { // loop for all particles in array
279 // Track to jet reordering
280 if(injet[jpart] == idx[kj]){
281 injetOk[jpart] = kj;
282 }
283 // Check if the particle belongs to the jet and add the ref
284 if(injetOk[jpart] == kj && fCalTrkEvent->GetCalTrkTrack(jpart)->GetCutFlag() == 1) {
285 jet.AddTrack(fCalTrkEvent->GetCalTrkTrack(jpart)->GetTrackObject());
286 }
287 }
288
289 AddJet(jet);
290
291 }
292
293 //delete
294 delete[] ptT;
295 delete[] etaT;
296 delete[] phiT;
297 delete[] injet;
298 delete[] injetOk;
299 delete[] areaJet;
300 delete[] areaJetOk;
301
302}
303
304//-----------------------------------------------------------------------
305void AliUA1JetFinder::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
306 Float_t* const etJet,Float_t* const etaJet, Float_t* const phiJet,
307 Float_t* const etallJet, Int_t* const ncellsJet)
308{
309 // Dump lego
310 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
311 const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory
312
313 const Int_t nBinEta = header->GetLegoNbinEta();
314 const Int_t nBinPhi = header->GetLegoNbinPhi();
315 if((nBinPhi*nBinEta)>nBinsMax){
316 AliError("Too many bins of the ETA-PHI histogram");
317 }
318
319 Float_t etCell[nBinsMax] = {0.}; //! Cell Energy
320 Float_t etaCell[nBinsMax] = {0.}; //! Cell eta
321 Float_t phiCell[nBinsMax] = {0.}; //! Cell phi
322 Short_t flagCell[nBinsMax] = {0}; //! Cell flag
323
324 Int_t nCell = 0;
325 TAxis* xaxis = fLego->GetXaxis();
326 TAxis* yaxis = fLego->GetYaxis();
327 Float_t e = 0.0;
328 for (Int_t i = 1; i <= nBinEta; i++) {
329 for (Int_t j = 1; j <= nBinPhi; j++) {
330 e = fLego->GetBinContent(i,j);
331 if (e < 0.0) continue; // don't include this cells
332 Float_t eta = xaxis->GetBinCenter(i);
333 Float_t phi = yaxis->GetBinCenter(j);
334 etCell[nCell] = e;
335 etaCell[nCell] = eta;
336 phiCell[nCell] = phi;
337 flagCell[nCell] = 0; //default
338 nCell++;
339 }
340 }
341 // Parameters from header
342 Float_t minmove = header->GetMinMove();
343 Float_t maxmove = header->GetMaxMove();
344 Float_t rc = header->GetRadius();
345 Float_t etseed = header->GetEtSeed();
346 // Tmp array of jets form algoritm
347 Float_t etaAlgoJet[kMaxJets] = {0.0};
348 Float_t phiAlgoJet[kMaxJets] = {0.0};
349 Float_t etAlgoJet[kMaxJets] = {0.0};
350 Int_t ncellsAlgoJet[kMaxJets] = {0};
351
352 // Run algorithm//
353
354 // Sort cells by et
355 Int_t index[nBinsMax];
356 TMath::Sort(nCell, etCell, index);
357 // variable used in centroide loop
358 Float_t eta = 0.0;
359 Float_t phi = 0.0;
360 Float_t eta0 = 0.0;
361 Float_t phi0 = 0.0;
362 Float_t etab = 0.0;
363 Float_t phib = 0.0;
364 Float_t etas = 0.0;
365 Float_t phis = 0.0;
366 Float_t ets = 0.0;
367 Float_t deta = 0.0;
368 Float_t dphi = 0.0;
369 Float_t dr = 0.0;
370 Float_t etsb = 0.0;
371 Float_t etasb = 0.0;
372 Float_t phisb = 0.0;
373 Float_t dphib = 0.0;
374
375 for(Int_t icell = 0; icell < nCell; icell++)
376 {
377 Int_t jcell = index[icell];
378 if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
379 if(flagCell[jcell] != 0) continue; // if cell was used before
380
381 eta = etaCell[jcell];
382 phi = phiCell[jcell];
383 eta0 = eta;
384 phi0 = phi;
385 etab = eta;
386 phib = phi;
387 ets = etCell[jcell];
388 etas = 0.0;
389 phis = 0.0;
390 etsb = ets;
391 etasb = 0.0;
392 phisb = 0.0;
393 for(Int_t kcell =0; kcell < nCell; kcell++)
394 {
395 Int_t lcell = index[kcell];
396 if(lcell == jcell) continue; // cell itself
397 if(flagCell[lcell] != 0) continue; // cell used before
398 if(etCell[lcell] > etCell[jcell]) continue; // can this happen
399 //calculate dr
400 deta = etaCell[lcell] - eta;
401 dphi = TMath::Abs(phiCell[lcell] - phi);
402 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
403 dr = TMath::Sqrt(deta * deta + dphi * dphi);
404 if(dr <= rc)
405 {
406 // calculate offset from initiate cell
407 deta = etaCell[lcell] - eta0;
408 dphi = phiCell[lcell] - phi0;
409 if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
410 if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
411 etas = etas + etCell[lcell]*deta;
412 phis = phis + etCell[lcell]*dphi;
413 ets = ets + etCell[lcell];
414 //new weighted eta and phi including this cell
415 eta = eta0 + etas/ets;
416 phi = phi0 + phis/ets;
417 // if cone does not move much, just go to next step
418 dphib = TMath::Abs(phi - phib);
419 if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
420 dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
421 if(dr <= minmove) break;
422 // cone should not move more than max_mov
423 dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
424 if(dr > maxmove){
425 eta = etab;
426 phi = phib;
427 ets = etsb;
428 etas = etasb;
429 phis = phisb;
430 } else { // store this loop information
431 etab=eta;
432 phib=phi;
433 etsb = ets;
434 etasb = etas;
435 phisb = phis;
436 }
437 } // inside cone
438 }//end of cells loop looking centroide
439
440 // Avoid cones overloap (to be implemented in the future)
441
442 // Flag cells in Rc, estimate total energy in cone
443 Float_t etCone = 0.0;
444 Int_t nCellIn = 0;
445 rc = header->GetRadius();
446
447 for(Int_t ncell =0; ncell < nCell; ncell++)
448 {
449 if(flagCell[ncell] != 0) continue; // cell used before
450 //calculate dr
451 deta = etaCell[ncell] - eta;
452 dphi = phiCell[ncell] - phi;
453 if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
454 if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
455 dr = TMath::Sqrt(deta * deta + dphi * dphi);
456 if(dr <= rc){ // cell in cone
457 flagCell[ncell] = -1;
458 etCone+=etCell[ncell];
459 nCellIn++;
460 }
461 }
462
463 // Select jets with et > background
464 // estimate max fluctuation of background in cone
465 Double_t ncellin = (Double_t)nCellIn;
466 Double_t ntcell = (Double_t)nCell;
467 Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
468 // min cone et
469 Double_t etcmin = etCone ; // could be used etCone - etmin !!
470 //decisions !! etbmax < etcmin
471
472 for(Int_t mcell =0; mcell < nCell; mcell++){
473 if(flagCell[mcell] == -1){
474 if(etbmax < etcmin)
475 flagCell[mcell] = 1; //flag cell as used
476 else
477 flagCell[mcell] = 0; // leave it free
478 }
479 }
480 //store tmp jet info !!!
481 if(etbmax < etcmin) {
482 if(nJets<kMaxJets){
483 etaAlgoJet[nJets] = eta;
484 phiAlgoJet[nJets] = phi;
485 etAlgoJet[nJets] = etCone;
486 ncellsAlgoJet[nJets] = nCellIn;
487 nJets++;
488 }
489 else{
490 AliError(Form("Too many jets (> %d) found by UA1JetFinder, adapt your cuts",kMaxJets));
491 break;
492 }
493 }
494 } // end of cells loop
495
496 //reorder jets by et in cone
497 //sort jets by energy
498 for(Int_t p = 0; p < nJets; p++)
499 {
500 etaJet[p] = etaAlgoJet[p];
501 phiJet[p] = phiAlgoJet[p];
502 etJet[p] = etAlgoJet[p];
503 etallJet[p] = etAlgoJet[p];
504 ncellsJet[p] = ncellsAlgoJet[p];
505 }
506
507}
508
509//-----------------------------------------------------------------------
510void AliUA1JetFinder::Reset()
511{
512 fLego->Reset();
513 AliJetFinder::Reset();
514
515}
516
517//-----------------------------------------------------------------------
518void AliUA1JetFinder::WriteJHeaderToFile() const
519{
520 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
521 header->Write();
522
523}
524
525//-----------------------------------------------------------------------
526void AliUA1JetFinder::Init()
527{
528
529 // initializes some variables
530 AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
531 // book lego
532 fLego = new TH2F("legoH","eta-phi",
533 header->GetLegoNbinEta(), header->GetLegoEtaMin(),
534 header->GetLegoEtaMax(), header->GetLegoNbinPhi(),
535 header->GetLegoPhiMin(), header->GetLegoPhiMax());
536
537}
538