]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - JETAN/AliCdfJetFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fix for coverity (AdC)
[u/mrichter/AliRoot.git] / JETAN / AliCdfJetFinder.cxx
... / ...
CommitLineData
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 usec, 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// Jet Finder based on CDF algortihm
20// Charged jet evolution and the underlying event in proton-antiproton collisions at 1.8 TeV
21// Physical Review D, vol. 65, Issue 9, id. 092002
22// http://www.phys.ufl.edu/~rfield/cdf/chgjet/chgjet_intro.html
23// Authors : Adrian.Sevcenco@cern.ch (adriansev@spacescience.ro )
24// Daniel.Felea@cern.ch (dfelea@spacescience.ro)
25// Ciprian.Mihai.Mitu@cern.ch (mcm@spacescience.ro)
26// magali.estienne@subatech.in2p3.fr &
27// alexandre.shabetai@cern.ch (Modification of the input object (reader/finder splitting))
28// ** 2011
29// Modified accordingly to reader/finder splitting and new handling of neutral information
30//---------------------------------------------------------------------
31
32#include <Riostream.h>
33#include <TMath.h>
34#include <TBits.h>
35#include <TFile.h>
36#include <TH1F.h>
37#include <TProfile.h>
38#include <TVector2.h>
39
40#include "AliAODJet.h"
41#include "AliJetFinder.h"
42#include "AliJetCalTrk.h"
43#include "AliCdfJetFinder.h"
44#include "AliCdfJetHeader.h"
45
46ClassImp(AliCdfJetFinder)
47
48//______________________________________________________________________________
49AliCdfJetFinder::AliCdfJetFinder():
50 AliJetFinder(),
51 fHistos(0),
52 fAODwrite(0),
53 fAODtracksWrite(0),
54 fAnalyseJets(0),
55 fNJets(0),
56 fNPart(0),
57 fNInC(0),
58 fNInN(0),
59 fRadius(0.7),
60 fMinJetParticles(1),
61 fJetPtCut(0.),
62 fVectParticle(NULL),
63 fVectJet(NULL),
64 fPtArray(NULL),
65 fIdxArray(NULL)
66{
67 // Default constructor
68}
69
70//______________________________________________________________________________
71AliCdfJetFinder::~AliCdfJetFinder()
72{
73 // destructor
74 Clean();
75
76}
77
78//______________________________________________________________________________
79void AliCdfJetFinder::CreateOutputObjects(TList * const histos)
80{
81 // Create the list of histograms. Only the list is owned.
82 fHistos = histos;
83
84 // gStyle->SetOptStat(11111111);
85
86 TH1F *h1 = new TH1F ("histo1", "Pt distribution of jets", 200, 0,200); // 1GeV/bin
87 h1->SetStats(kTRUE);
88 h1->GetXaxis()->SetTitle("P_{T} of jets");
89 h1->GetYaxis()->SetTitle("Number of jets");
90 h1->GetXaxis()->SetTitleColor(1);
91 h1->SetMarkerStyle(kFullCircle);
92 fHistos->Add(h1);
93
94 TH1F *h2 = new TH1F ("histo2", "Eta distribution of jets", 240, -1.2,1.2); // 1 unit of rapidity / 100 bin
95 h2->SetStats(kTRUE);
96 h2->GetXaxis()->SetTitle("Eta of jets");
97 h2->GetYaxis()->SetTitle("Number of jets");
98 h2->GetXaxis()->SetTitleColor(1);
99 h2->SetMarkerStyle(kFullCircle);
100 fHistos->Add(h2);
101
102 TH1F *h3 = new TH1F ("histo3", "Phi distribution of jets", 400, -4,4);
103 h3->SetStats(kTRUE);
104 h3->GetXaxis()->SetTitle("Phi of jets");
105 h3->GetYaxis()->SetTitle("Number of jets");
106 h3->GetXaxis()->SetTitleColor(1);
107 h3->SetMarkerStyle(kFullCircle);
108 fHistos->Add(h3);
109
110 TH1F *h4 = new TH1F ("histo4", "Multiplicity of jets", 40, 0,40); // 1 unit of multiplicity /bin
111 h4->SetStats(kTRUE);
112 h4->GetXaxis()->SetTitle("Particles in jets");
113 h4->GetYaxis()->SetTitle("Number of jets");
114 h4->GetXaxis()->SetTitleColor(1);
115 h4->SetMarkerStyle(kFullCircle);
116 fHistos->Add(h4);
117
118 TH1F *h5 = new TH1F ("histo5", "Distribution of jets in events", 100, 0,100);
119 h5->SetStats(kTRUE);
120 h5->GetXaxis()->SetTitle("Number of jets");
121 h5->GetYaxis()->SetTitle("Number of events");
122 h5->GetXaxis()->SetTitleColor(1);
123 h5->SetMarkerStyle(kFullCircle);
124 fHistos->Add(h5);
125
126 TH1F *h6 = new TH1F ("histo6", "Jet1 Charged Multiplicity Distribution", 30, 0,30);
127 h6->SetStats(kTRUE);
128 h6->GetXaxis()->SetTitle("N_{chg}");
129 h6->GetYaxis()->SetTitle("Number of jets");
130 h6->GetXaxis()->SetTitleColor(1);
131 h6->SetMarkerStyle(kFullCircle);
132 fHistos->Add(h6);
133
134 TProfile * h7 = new TProfile ("histo7","N_{chg}(jet1) vs P_{T}(charged jet1)", 200, 0. ,200. , 0.,200. ) ;
135 h7->SetStats(kTRUE);
136 h7->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
137 h7->GetYaxis()->SetTitle("<N_{chg}(jet1)> in 1 GeV/c bin");
138 h7->GetXaxis()->SetTitleColor(1);
139 h7->SetMarkerStyle(kFullCircle);
140 fHistos->Add(h7);
141
142 TH1F *h8 = new TH1F ("histo8", "Charge momentum distribution for leading jet", 120, 0 , 1.2);
143 h8->SetStats(kTRUE);
144 h8->GetXaxis()->SetTitle("Jets");
145 h8->GetYaxis()->SetTitle("Particle distribution");
146 h8->GetXaxis()->SetTitleColor(1);
147 h8->SetMarkerStyle(kFullCircle);
148 fHistos->Add(h8);
149
150 TProfile *h9 = new TProfile ("histo9", "N_{chg} vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 20 );
151 h9->SetStats(kTRUE);
152 h9->GetXaxis()->SetTitle("#Delta#phi (degrees)");
153 h9->GetYaxis()->SetTitle("<N_{chg}> in 3.6 degree bin");
154 h9->GetXaxis()->SetTitleColor(1);
155 h9->SetMarkerStyle(kFullCircle);
156 fHistos->Add(h9);
157
158 TProfile *h10 = new TProfile ("histo10", "P_{T} sum vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 100 );
159 h10->SetStats(kTRUE);
160 h10->GetXaxis()->SetTitle("#Delta#phi (degrees)");
161 h10->GetYaxis()->SetTitle("<P_{T} sum> in 3.6 degree bin");
162 h10->GetXaxis()->SetTitleColor(1);
163 h10->SetMarkerStyle(kFullCircle);
164 fHistos->Add(h10);
165
166 TH1F *h11 = new TH1F ("histo11", " \"Transverse\" Pt Distribution ", 70, 0 , 14);
167 h11->SetStats(kTRUE);
168 h11->GetXaxis()->SetTitle("P_{T} (GeV/c)");
169 h11->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
170 h11->GetXaxis()->SetTitleColor(1);
171 h11->SetMarkerStyle(kFullCircle);
172 fHistos->Add(h11);
173
174 TH1F *h20 = new TH1F ("histo20", "Distribution of R in leading jet", 400, 0.,4.);
175 h20->SetStats(kTRUE);
176 h20->GetXaxis()->SetTitle("R [formula]");
177 h20->GetYaxis()->SetTitle("dN/dR");
178 h20->GetXaxis()->SetTitleColor(1);
179 h20->SetMarkerStyle(kFullCircle);
180 fHistos->Add(h20);
181
182 TProfile * h21 = new TProfile ("histo21","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 30. ) ;
183 h21->SetStats(kTRUE);
184 h21->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
185 h21->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
186 h21->GetXaxis()->SetTitleColor(1);
187 h21->SetMarkerStyle(kFullCircle);
188 fHistos->Add(h21);
189
190 TProfile * h22 = new TProfile ("histo22","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
191 h22->SetStats(kTRUE);
192 h22->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
193 h22->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
194 h22->GetXaxis()->SetTitleColor(1);
195 h22->SetMarkerStyle(kFullCircle);
196 fHistos->Add(h22);
197
198 TProfile * h21Toward = new TProfile ("histo21_toward","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
199 h21Toward->SetStats(kTRUE);
200 h21Toward->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
201 h21Toward->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
202 h21Toward->GetXaxis()->SetTitleColor(1);
203 h21Toward->SetMarkerStyle(kFullCircle);
204 fHistos->Add(h21Toward);
205
206 TProfile * h21Transverse = new TProfile ("histo21_transverse","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
207 h21Transverse->SetStats(kTRUE);
208 h21Transverse->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
209 h21Transverse->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
210 h21Transverse->GetXaxis()->SetTitleColor(1);
211 h21Transverse->SetMarkerStyle(kFullCircle);
212 fHistos->Add(h21Transverse);
213
214 TProfile * h21Away = new TProfile ("histo21_away","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
215 h21Away->SetStats(kTRUE);
216 h21Away->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
217 h21Away->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
218 h21Away->GetXaxis()->SetTitleColor(1);
219 h21Away->SetMarkerStyle(kFullCircle);
220 fHistos->Add(h21Away);
221
222 TProfile * h22Toward = new TProfile ("histo22_toward","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
223 h22Toward->SetStats(kTRUE);
224 h22Toward->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
225 h22Toward->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
226 h22Toward->GetXaxis()->SetTitleColor(1);
227 h22Toward->SetMarkerStyle(kFullCircle);
228 fHistos->Add(h22Toward);
229
230 TProfile * h22Transverse = new TProfile ("histo22_transverse","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
231 h22Transverse->SetStats(kTRUE);
232 h22Transverse->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
233 h22Transverse->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
234 h22Transverse->GetXaxis()->SetTitleColor(1);
235 h22Transverse->SetMarkerStyle(kFullCircle);
236 fHistos->Add(h22Transverse);
237
238 TProfile * h22Away = new TProfile ("histo22_away","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
239 h22Away->SetStats(kTRUE);
240 h22Away->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
241 h22Away->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
242 h22Away->GetXaxis()->SetTitleColor(1);
243 h22Away->SetMarkerStyle(kFullCircle);
244 fHistos->Add(h22Away);
245
246 TH1F *h23Toward = new TH1F ("histo23_toward","'Toward' Pt Distribution of charged particles", 200, 0., 14.);
247 h23Toward->SetStats(kTRUE);
248 h23Toward->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
249 h23Toward->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
250 h23Toward->GetXaxis()->SetTitleColor(1);
251 h23Toward->SetMarkerStyle(kFullCircle);
252 fHistos->Add(h23Toward);
253
254 TH1F *h23Transverse = new TH1F ("histo23_transverse","'Transverse' Pt Distribution of charged particles", 200, 0., 14.);
255 h23Transverse->SetStats(kTRUE);
256 h23Transverse->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
257 h23Transverse->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
258 h23Transverse->GetXaxis()->SetTitleColor(1);
259 h23Transverse->SetMarkerStyle(kFullCircle);
260 fHistos->Add(h23Transverse);
261
262 TH1F *h23Away = new TH1F ("histo23_away","'Away' Pt Distribution of charged particles", 200, 0., 14.);
263 h23Away->SetStats(kTRUE);
264 h23Away->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
265 h23Away->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
266 h23Away->GetXaxis()->SetTitleColor(1);
267 h23Away->SetMarkerStyle(kFullCircle);
268 fHistos->Add(h23Away);
269
270 TProfile * h24 = new TProfile ("histo24","Jet1 Size vs P_{T}(charged jet1)", 200, 0., 50. , 0., 0.5) ;
271 h24->SetStats(kTRUE);
272 h24->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
273 h24->GetYaxis()->SetTitle("<R(chgjet1)> in 1 GeV/c bin");
274 h24->GetXaxis()->SetTitleColor(1);
275 h24->SetMarkerStyle(kFullCircle);
276 fHistos->Add(h24);
277
278 TProfile * h25 = new TProfile ("histo25","Jet1 Size vs P_{T}(charged jet1)", 200, 0., 50. , 0., 0.5) ;
279 h25->SetStats(kTRUE);
280 h25->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
281 h25->GetYaxis()->SetTitle("<R(chgjet1)> in 1 GeV/c bin");
282 h25->GetXaxis()->SetTitleColor(1);
283 h25->SetMarkerStyle(kFullCircle);
284 fHistos->Add(h25);
285
286 TProfile *h26 = new TProfile ("histo26", "N_{chg} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0., 0.8);
287 h26->SetStats(kTRUE);
288 h26->GetXaxis()->SetTitle("Distance R");
289 h26->GetYaxis()->SetTitle("<N_{chg}> in 0.02 bin");
290 h26->GetXaxis()->SetTitleColor(1);
291 h26->SetMarkerStyle(kFullCircle);
292 fHistos->Add(h26);
293
294 TProfile *h27 = new TProfile ("histo27", "N_{chg} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0., 0.8);
295 h27->SetStats(kTRUE);
296 h27->GetXaxis()->SetTitle("Distance R");
297 h27->GetYaxis()->SetTitle("<N_{chg}> in 0.02 bin");
298 h27->GetXaxis()->SetTitleColor(1);
299 h27->SetMarkerStyle(kFullCircle);
300 fHistos->Add(h27);
301
302 TProfile *h28 = new TProfile ("histo28", "PT_{sum} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0.01, 10.);
303 h28->SetStats(kTRUE);
304 h28->GetXaxis()->SetTitle("Distance R");
305 h28->GetYaxis()->SetTitle("<PT_{sum} (GeV/c)> in 0.02 bin");
306 h28->GetXaxis()->SetTitleColor(1);
307 h28->SetMarkerStyle(kFullCircle);
308 fHistos->Add(h28);
309
310 TProfile *h29 = new TProfile ("histo29", "PT_{sum} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0.01, 10.);
311 h29->SetStats(kTRUE);
312 h29->GetXaxis()->SetTitle("Distance R");
313 h29->GetYaxis()->SetTitle("<PT_{sum} (GeV/c)> in 0.02 bin");
314 h29->GetXaxis()->SetTitleColor(1);
315 h29->SetMarkerStyle(kFullCircle);
316 fHistos->Add(h29);
317
318}
319
320//______________________________________________________________________________
321void AliCdfJetFinder::FindJets()
322{
323 // Jet Algorithm:
324 // * Order all charged particles according to their PT.
325 // * Start with the highest PT particle and include in the "jet" all particles within the "radius" R = 0.7
326 // (considering each particle in the order of decreasing PT and recalculating the centroid of the jet after
327 // each new particle is added to the jet).
328 // * Go to the next highest PT particle (not already included in a jet) and include in the "jet" all particles
329 // (not already included in a jet) within the radius R =0.7.
330 // * Continue until all particles are in a "jet".
331 if (fDebug) { printf("AliCDJetfinder::FindJets() %d \n", __LINE__ ); }
332 AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;
333
334 if (header)
335 {
336 fDebug = header->GetDebug();
337 fAODwrite = header->IsAODwrite() ; // write jets to AOD
338 fAODtracksWrite = header->IsAODtracksWrite() ; // write jet tracks to AOD
339 fRadius = header->GetRadius(); // get Radius from jet finder header
340 fMinJetParticles = header->GetMinPartJet (); // get minimum multiplicity of an jet
341 fJetPtCut = header->GetJetPtCut (); // get minimum of jet pt
342 fAnalyseJets = header->GetAnalyseJets(); // get analyse jet
343 }
344 else
345 { cout << "Header not found" << endl; return; }
346
347 InitData();
348
349 if (!fNPart) {
350 if (fDebug) {
351 cout << "entries = 0 ; Event empty !!!" << endl ;
352 }
353 // no need to call clean, InitData does not
354 // create pointers if npart == 0
355 return;
356 } // if event empty then exit
357
358 FindCones();
359
360 ComputeConesWeight();
361
362 if (fAODwrite) {
363 if(fDebug)cout << "Writing AOD" << endl ;
364 WriteJets();
365 }
366
367 if (fAnalyseJets) AnalizeJets();
368
369 Clean();
370
371}
372
373//______________________________________________________________________________
374void AliCdfJetFinder::InitData()
375{
376 // initialisation of variables and data members
377
378 if( fHeader->GetDebug() && fCalTrkEvent->GetNCalTrkTracks() == 0) { cout << "No charged tracks found" << endl; }
379
380 fNPart = fCalTrkEvent->GetNCalTrkTracks() ;
381
382 if ( fCalTrkEvent->GetNCalTrkTracks() ) { return; } // if event empty then exit
383
384 fVectParticle = new varContainer* [fNPart]; // container for Particles
385
386 fPtArray = new Double_t [fCalTrkEvent->GetNCalTrkTracks()] ;
387 fIdxArray = new Int_t [fCalTrkEvent->GetNCalTrkTracks()] ; // index array of sorted pts
388
389 // initialisation of momentum and index arrays
390 for ( Int_t i = 0 ; i <fCalTrkEvent->GetNCalTrkTracks() ; i++ )
391 {// SORTING STEP :: fPtArray with data from CalTrkTracks
392
393 // INITIALISATION of local arrays for temporary storage
394 varContainer *aParticle = new varContainer;
395 aParticle->pt = fCalTrkEvent->GetCalTrkTrack(i)->GetPt();
396 aParticle->eta = fCalTrkEvent->GetCalTrkTrack(i)->GetEta();
397 aParticle->phi = TVector2::Phi_mpi_pi ( fCalTrkEvent->GetCalTrkTrack(i)->GetPhi() ); // normalize to -pi,pi
398 aParticle->njet = -999;
399
400 fVectParticle[i] = aParticle; // vector of Particles
401
402 // initializing arrays
403 fIdxArray [i] = -999 ;
404 fPtArray [i] = aParticle->pt ;
405 }
406
407 TMath::Sort ( fNPart, fPtArray, fIdxArray ) ; // get a sorted array of indexes
408
409}
410
411//______________________________________________________________________________
412void AliCdfJetFinder::FindCones()
413{
414 // parsing of particles in event and estlabish jets (label them with jet index)
415
416 Double_t ptSeed = 0. , etaSeed = 0. , phiSeed = 0. ; // leading particle params
417 Double_t pttmp = 0. , etatmp = 0. , phitmp = 0. ; // temporary variables to be used in various calculations
418 Double_t deta = 0. , dphi = 0. , dcomputed = 0. ;
419 Bool_t injet = 0 ;
420
421 fNJets = -1 ; // n jets in this event
422 Int_t idxPtSort = -1 ; // index of array of sorted pt indexes
423
424 if (fDebug) { cout << "\n\n\n\n\n\n------------------\nBegin Event Analysis\n------------------\n\n" << endl ;}
425
426 if(fDebug)cout << "fNPart = " << fNPart << endl;
427
428 TBits lkupTable ( fNPart ) ; // bit container ; 1-to-1 corespondence with fIdxArray
429
430 while ( lkupTable.CountBits() != (UInt_t)fNPart )
431 { // loop over particles in event until all flags are set
432 UInt_t firstnonflagged = lkupTable.FirstNullBit() ; // set the index to the first NON flagged bit ; less conditions
433
434 if(fDebug)cout << "\n\nfirst_non_flagged : " << firstnonflagged << endl;
435
436 ++fNJets; // incrementing the jet counter
437 if (fDebug) { printf("JET %d \n", fNJets); }
438
439 ptSeed = 0. ; etaSeed = 0. ; phiSeed = 0. ; // reseting leading particle params
440
441 for ( UInt_t ipart = firstnonflagged ; ipart < (UInt_t)fNPart ; ipart++ )
442 {// iteration over particles in event
443 // the loop is done over sorted array of pt
444 idxPtSort = fIdxArray[ipart] ; // index of particle ! fIdxArray is an index list pt sorted
445
446 if ( lkupTable.TestBitNumber(ipart) ) { continue; } // if 4vector is already flagged skip it
447
448 //init computed and used vars
449 pttmp = 0. ; etatmp = 0. ; phitmp = 0. ;
450 deta = 0. ; dphi = 0. ; dcomputed = 0. ; injet = 0 ;
451
452 //taking info from fVectParticle ;
453 pttmp = fVectParticle[idxPtSort]->pt ;
454 etatmp = fVectParticle[idxPtSort]->eta ;
455 phitmp = fVectParticle[idxPtSort]->phi ;
456
457 if ( ipart == firstnonflagged )
458 {// this is first particle in event; leading particle
459 // begin the search around this particle in a fRadius
460
461 // CENTRE OF THE JET
462 ptSeed = pttmp ; etaSeed = etatmp ; phiSeed = phitmp ; // seeding the jet with first particle idxPtSort
463
464 lkupTable.SetBitNumber ( ipart ) ; // flag the index of particle in lkup_table
465 fVectParticle[idxPtSort]->njet = fNJets ; // associate particle with current jet number
466
467 if (fDebug) { printf("\nLeading particle :: particle index = %d ; at sorted index %d ; in jet %d \n", idxPtSort, ipart, fNJets); }
468 if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pttmp, etatmp, phitmp) ; }
469 if (fDebug) { lkupTable.Print() ;}
470
471 continue ; // skip to next particle
472 }
473
474 // condition to be in jet
475 deta = etatmp - etaSeed ;
476 dphi = TVector2::Phi_mpi_pi (phitmp - phiSeed) ; // computing dphi and normalizing to (0,2pi) interval in one step
477
478 dcomputed = TMath::Hypot(deta, dphi) ; // Distance(fRadius) to (eta,phi) seed
479
480 injet = ( ( fRadius - dcomputed ) >= 0.000000001 ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0
481
482 if ( injet )
483 { // calculus of jet variables
484 lkupTable.SetBitNumber ( ipart ) ; // flag the index of particle in lkup_table
485 fVectParticle[idxPtSort]->njet = fNJets ; // setting in particle list the associated jet
486
487 if (fDebug) { printf("\njet particle :: particle index = %d ; at sorted index %d ; in jet %d ; found at radius %g ; \n", idxPtSort, ipart, fNJets, dcomputed); }
488 if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pttmp, etatmp, phitmp) ; }
489 if (fDebug) { lkupTable.Print() ;}
490
491 continue ; // skip to next particle
492 }
493
494 }
495 // end of iteration over event; one jet definition of content ; jet parameters to be computed later
496 }
497
498}
499
500//______________________________________________________________________________
501void AliCdfJetFinder::ComputeConesWeight()
502{
503 // computing of jets Pt, Eta and Phi (centre of weight in (eta,phi) plane)
504 // rescan the vector of particles by identify them by asociate jet number for computing of weight centre
505
506 // JET CONTAINER
507 fVectJet = new varContainer* [fNJets]; // container for Jets
508
509 Double_t ptJet, ptJet2 , etaJet , phiJet ; Int_t npartJet ;
510 Double_t pttmp = 0. , etatmp = 0. , phitmp = 0. ; // temporary variables to be used in various calculations
511 Int_t idxPtSort = -999 ; // index of array of sorted pt indexes
512
513 for( Int_t jet = 0 ; jet < fNJets ; jet++ )
514 {
515 if (fDebug) { printf("\n\n--- Computing weight of Jet %d \n", jet ); }
516 npartJet = 0 ; ptJet = 0. ; etaJet = 0. ; phiJet = 0. ; // reset variables for a new computation
517
518 for ( Int_t ipart = 0 ; ipart < fNPart ; ipart++ )
519 {// iteration over particles in event
520 // the loop is done over sorted array of pt
521 idxPtSort = fIdxArray[ipart] ; // index of particle ! fIdxArray is an index list pt sorted
522
523 if ( fVectParticle[idxPtSort]->njet == jet )
524 {
525 ++npartJet; // incrementing the counter of jet particles
526
527 //taking info from fVectParticle ;
528 pttmp = fVectParticle[idxPtSort]->pt ;
529 etatmp = fVectParticle[idxPtSort]->eta ;
530 phitmp = TVector2::Phi_mpi_pi (fVectParticle[idxPtSort]->phi) ;
531
532 // jet_new_angular_coordinate = jet_old_angular_coordinate * jet_old_pt / jet_new_pt +
533 // part[i]_angular_coordinate * part[i]_pt/jet_new_pt
534
535 ptJet2 = ptJet + pttmp ;
536
537 etaJet = etaJet * ptJet / ptJet2 + etatmp * pttmp / ptJet2 ;
538 phiJet = phiJet * ptJet / ptJet2 + phitmp * pttmp / ptJet2 ;
539
540 ptJet = ptJet2 ;
541
542 }
543 // add a particle and recalculation of centroid
544 }
545 // end of 1 jet computation
546
547 varContainer *aJet = new varContainer; // Jet container
548 aJet->pt = ptJet; aJet->eta = etaJet; aJet->phi = phiJet; aJet->njet = npartJet; // setting jet vars in container
549 fVectJet[jet] = aJet; // store the number of the jet(fNJets) and increment afterwards
550
551 if (fDebug) { printf ("=== current jet %d : npartjet= %d ; pt_jet= %g ; eta_jet = %g ; phi_jet = %g \n\n\n",
552 jet, npartJet, ptJet, etaJet, phiJet ) ; }
553
554 }
555 //end loop over jets
556
557}
558
559//______________________________________________________________________________
560void AliCdfJetFinder::WriteJets()
561{
562 // Writing AOD jets and AOD tracks
563
564 for( Int_t jetnr = 0 ; jetnr < fNJets ; jetnr++ )
565 {
566 Double_t pt = 0., eta = 0., phi = 0., // jet variables
567 px = 0., py = 0., pz = 0., en = 0.; // convert to 4-vector
568 pt = fVectJet[ jetnr ]->pt ; // pt of jet
569 eta = fVectJet[ jetnr ]->eta ; // eta of jet
570 phi = fVectJet[ jetnr ]->phi ; // phi of jet
571
572 px = pt * TMath::Cos ( phi ) ;
573 py = pt * TMath::Sin ( phi ) ;
574 pz = pt / TMath::Tan ( 2.0 * TMath::ATan ( TMath::Exp ( -eta ) ) ) ;
575 en = TMath::Sqrt ( px * px + py * py + pz * pz );
576
577 AliAODJet jet (px, py, pz, en);
578
579
580 if (fDebug) jet.Print("");
581
582 if (fAODtracksWrite)
583 {
584 for ( Int_t jetTrack = 0; jetTrack < fCalTrkEvent->GetNCalTrkTracks(); jetTrack++ )
585 {
586 // The first if condition below has to be checked
587 if ( fVectParticle[jetTrack]->njet == jetnr ) { jet.AddTrack(fCalTrkEvent->GetCalTrkTrack(jetTrack)->GetTrackObject()) ; }
588 }
589 }
590 // tracks REFs written in AOD
591 AddJet(jet);
592
593 }
594 //jets vector parsed and written to AOD
595}
596
597//______________________________________________________________________________
598void AliCdfJetFinder::AnalizeJets()
599{
600 // analyzing of jets and filling of histograms
601
602 const Double_t kPI = TMath::Pi();
603
604 //persistent pointer to histo20
605 TH1F *hR = (TH1F*)fHistos->FindObject("histo20");
606
607 Int_t *jetsptidx = 0; // sorted array of jets pt
608 Double_t *jetspt = 0; // array of jets pts
609 Int_t leadingjetindex = -1 ; // index of leading jet from fVectJet
610 Int_t partleadjet = 0 ; // number of particles in leading jet
611 Double_t ptleadjet = 0. ; // pt of leading jet
612 Double_t etaleadjet = 0. ; // eta of leading jet
613 Double_t phileadjet = 0. ; // phi of leading jet
614
615 jetsptidx = new Int_t [fNJets] ;
616 jetspt = new Double_t [fNJets] ;
617
618 //________________________________________________________________________________________
619 // Jet sorting and finding the leading jet that coresponds to cuts in pt and multiplicity
620 //________________________________________________________________________________________
621
622 // filing the idx_ptjets array
623 if (fDebug) printf("List of unsorted jets:\n");
624 for( Int_t i = 0 ; i < fNJets ; i++ )
625 {
626 jetsptidx [i] = 0 ;
627 jetspt [i] = fVectJet[i]->pt ;
628 if (fDebug) { cout << " jet found: " << i << " npartjet=" << fVectJet[i]->njet << " ; jets_pt = " << jetspt[i] << endl; }
629 }
630
631 TMath::Sort ( fNJets, jetspt , jetsptidx ) ; // sorting pt of jets
632
633 // selection of leading jet
634 // looping over jets searching for __first__ one that coresponds to cuts
635 for( Int_t i = 0 ; i < fNJets ; i++ )
636 {
637 if ( ( fVectJet[ jetsptidx[i] ]->njet >= fMinJetParticles ) && ( fVectJet[ jetsptidx[i] ]->pt >= fJetPtCut ) )
638 {
639 leadingjetindex = jetsptidx[i] ;
640 partleadjet = fVectJet[ leadingjetindex ]->njet ; // number of particles in leading jet
641 ptleadjet = fVectJet[ leadingjetindex ]->pt ; // pt of leading jet
642 etaleadjet = fVectJet[ leadingjetindex ]->eta ; // eta of leading jet
643 phileadjet = fVectJet[ leadingjetindex ]->phi ; // phi of leading jet
644
645 if (fDebug)
646 { printf("Leading jet %d : npart= %d ; pt= %g ; eta = %g ; phi = %g \n", leadingjetindex, partleadjet, ptleadjet, etaleadjet, phileadjet ); }
647
648 break ;
649 }
650 }
651 // end of selection of leading jet
652
653
654
655 //////////////////////////////////////////////////
656 //// Computing of values used in histograms
657 //////////////////////////////////////////////////
658
659 //___________________________________________________________________________
660 // pt_sum of all particles in event
661 //___________________________________________________________________________
662 if (fDebug) cout << "Computing sum of pt in event" << endl ;
663 Double_t ptsumevent = 0.;
664 for ( Int_t i = 0 ; i< fNPart ; i++ ) { ptsumevent += fVectParticle[i]->pt ; }
665 if (fDebug) printf ("Sum of all Pt in event : pt_sum_event = %g", ptsumevent) ;
666
667 //___________________________________________________________________________
668 // Filling an array with indexes of leading jet particles
669 //___________________________________________________________________________
670 Int_t * idxpartLJ = new Int_t [partleadjet] ;
671 Int_t counterpartleadjet = 0;
672
673 if (fDebug) cout << "Filling an array with indexes of leading jet particles" << endl;
674
675 for( Int_t i = 0 ; i < fNPart ; i++ )
676 {
677 if ( fVectParticle[i]->njet == leadingjetindex )
678 { idxpartLJ[counterpartleadjet++] = i ; }
679 }
680
681 if ( (counterpartleadjet-1) > partleadjet ) { cout << " Counter_part_lead_jet > part_leadjet !!!!" << endl;}
682
683
684 //___________________________________________________________________________
685 // Calculus of part distribution in leading jet
686 //___________________________________________________________________________
687 Double_t z = 0. ;
688 Double_t *zpartljet = new Double_t [ partleadjet ] ; // array of z of particles in leading jet
689
690 if (fDebug) cout << "Entering loop of calculus of part distribution in leading jet" << endl ;
691
692 for( Int_t j = 0 ; j < partleadjet ; j++ )
693 {
694 Double_t zj = fVectParticle[idxpartLJ[j]]->pt ;
695 z = zj / ptleadjet ;
696 zpartljet [j] = z ;
697 if (fDebug) cout << "idx_leadjet_part[j] = " << idxpartLJ[j]
698 << " p of particle = " << zj
699 << " pt lead jet = " << ptleadjet
700 << " Z = " << z << endl;
701 }
702
703
704 //___________________________________________________________________________
705 // array of delta phi's between phi of particles and leading jet phi
706 //___________________________________________________________________________
707 if (fDebug) cout << "array of delta phi's between phi of particles and leading jet phi" << endl;
708 Double_t dphipartLJ = 0. ;
709 Double_t *dphipartljet = new Double_t [fNPart];
710 for( Int_t part = 0 ; part < fNPart ; part++ )
711 {
712 dphipartLJ = fVectParticle[part]->phi - phileadjet ;
713 dphipartLJ = TVector2::Phi_mpi_pi (dphipartLJ) ; // restrict the delta phi to (-pi,pi) interval
714 dphipartljet [part] = dphipartLJ ;
715 if (fDebug) printf("part= %d ; dphi_partLJ = %g \n", part, dphipartLJ );
716 }
717
718
719 //______________________________________________________________________________
720 // Pt distribution for all particles
721 //______________________________________________________________________________
722 TH1F * hpt = (TH1F*)fHistos->FindObject("histo11");
723 if ( hpt ) { for ( Int_t i = 0 ; i < fNPart ; i++ ) { hpt->Fill( fVectParticle[i]->pt ); } }
724
725 //___________________________________________________________________________
726 // Recomputing of radius of particles in leading jet
727 //___________________________________________________________________________
728 if (fDebug) { printf(" Searching particles with jet index %d\n", leadingjetindex); }
729
730 Double_t ddeta = 0. , ddphi = 0. , rpart = 0. ;
731
732 for( Int_t j = 0 ; j < partleadjet ; j++ )
733 {
734 ddeta = etaleadjet - fVectParticle[idxpartLJ[j]]->eta;
735
736 Double_t phitmp = fVectParticle[idxpartLJ[j]]->phi ;
737
738 ddphi = TVector2::Phi_mpi_pi ( phileadjet - phitmp ) ; // restrict the delta phi to (-pi,pi) interval
739
740 rpart = TMath::Hypot (ddeta, ddphi) ;
741
742 if (fDebug) printf ("Particle %d with Re-Computed radius = %f ", idxpartLJ[j], rpart) ;
743 if ( (rpart - fRadius) >= 0.00000001 )
744 { if (fDebug) printf (" bigger than selected radius of %f\n", fRadius ); }
745 else
746 { if (fDebug) printf ("\n") ; }
747
748 if (hR) hR->Fill(rpart);
749
750 }
751
752
753
754 //_______________________________________________________________________
755 // Computing of radius that contain 80% of Leading Jet ( PT and multiplicity )
756 //_______________________________________________________________________
757 Double_t corepartleadjet = 0.8 * partleadjet ;
758 Double_t coreptleadjet = 0.8 * ptleadjet ;
759 Int_t countercorepart = 0 ;
760 Double_t countercorept = 0. ;
761 Int_t sortedindex = -1 ;
762
763 TProfile * hprof24 = (TProfile*)fHistos->FindObject("histo24");
764 TProfile * hprof25 = (TProfile*)fHistos->FindObject("histo25");
765 TProfile * hprof26 = (TProfile*)fHistos->FindObject("histo26");
766 TProfile * hprof27 = (TProfile*)fHistos->FindObject("histo27");
767 TProfile * hprof28 = (TProfile*)fHistos->FindObject("histo28");
768 TProfile * hprof29 = (TProfile*)fHistos->FindObject("histo29");
769
770
771 if ((hprof24) && (hprof25) && (hprof26) && (hprof27) && (hprof28) && (hprof29) )
772 {
773 for( Int_t part = 0 ; part < fNPart ; part++ )
774 {
775 Double_t pttmp = 0. ; Double_t etatmp = 0. ; Double_t phitmp = 0. ; // temporary variables
776 Double_t dpart = 0. ;
777 sortedindex = fIdxArray[part] ;
778
779 if ( fVectParticle [ sortedindex ]->njet == leadingjetindex )
780 {
781 pttmp = fVectParticle[sortedindex]->pt ;
782 etatmp = fVectParticle[sortedindex]->eta ;
783 phitmp = fVectParticle[sortedindex]->phi ;
784
785 ++countercorepart ;
786 countercorept += pttmp ;
787
788 dpart = TMath::Hypot ( etaleadjet - etatmp, TVector2::Phi_mpi_pi (phileadjet - phitmp) ) ;
789
790 if ( countercorepart <= corepartleadjet ) { hprof24->Fill(ptleadjet, dpart); }
791 if ( countercorept <= coreptleadjet ) { hprof25->Fill(ptleadjet, dpart); }
792
793 if (ptleadjet > 5.) { hprof26->Fill(dpart, countercorepart); hprof28->Fill(dpart, countercorept); }
794 if (ptleadjet > 30.) { hprof27->Fill(dpart, countercorepart); hprof29->Fill(dpart, countercorept); }
795
796 }
797 }
798 }
799
800 TH1F *hjetpt = (TH1F*)fHistos->FindObject("histo1");
801 TH1F *hjeteta = (TH1F*)fHistos->FindObject("histo2");
802 TH1F *hjetphi = (TH1F*)fHistos->FindObject("histo3");
803 TH1F *hjetnjet = (TH1F*)fHistos->FindObject("histo4");
804
805 for( Int_t jet = 0 ; jet < fNJets ; jet++ )
806 {
807 if (hjetpt) hjetpt ->Fill ( fVectJet[jet]->pt ) ;
808 if (hjeteta) hjeteta ->Fill ( fVectJet[jet]->eta ) ;
809 if (hjetphi) hjetphi ->Fill ( fVectJet[jet]->phi ) ;
810 if (hjetnjet) hjetnjet ->Fill ( fVectJet[jet]->njet ) ;
811 }
812
813 TH1F *hjets = (TH1F*)fHistos->FindObject("histo5");
814 if (hjets) hjets->Fill(fNJets);
815
816 TH1F *hleadpart = (TH1F*)fHistos->FindObject("histo6");
817 if (hleadpart) hleadpart->Fill(partleadjet);
818
819 TProfile * hprof = (TProfile*)fHistos->FindObject("histo7");
820 if (hprof) hprof->Fill(ptleadjet,partleadjet);
821
822 TH1F *hMD = (TH1F*)fHistos->FindObject("histo8");
823 for( Int_t k = 0 ; k < partleadjet ; k++)
824 { hMD->Fill( zpartljet[k] ); }
825
826 TProfile * hphi = (TProfile*)fHistos->FindObject("histo9");
827 for( Int_t k = 0 ; k < partleadjet ; k++)
828 { hphi->Fill( TMath::RadToDeg() * dphipartljet [k] , fNPart ) ; }
829
830 TProfile * htpd = (TProfile*)fHistos->FindObject("histo10");
831 for( Int_t k = 0 ; k < fNPart ; k++)
832 { htpd->Fill( TMath::RadToDeg() * dphipartljet [k] , ptsumevent ) ; }
833
834
835 TProfile * hprof1 = (TProfile*)fHistos->FindObject("histo21");
836 if (hprof1) hprof1->Fill(ptleadjet, fNPart);
837
838 TProfile * hprof2 = (TProfile*)fHistos->FindObject("histo22");
839 if (hprof2) hprof2->Fill(ptleadjet, ptsumevent);
840
841 TProfile * hprof1toward = (TProfile*)fHistos->FindObject("histo21_toward");
842 TProfile * hprof1transverse = (TProfile*)fHistos->FindObject("histo21_transverse");
843 TProfile * hprof1away = (TProfile*)fHistos->FindObject("histo21_away");
844 TProfile * hprof2toward = (TProfile*)fHistos->FindObject("histo22_toward");
845 TProfile * hprof2transverse = (TProfile*)fHistos->FindObject("histo22_transverse");
846 TProfile * hprof2away = (TProfile*)fHistos->FindObject("histo22_away");
847 TH1F * hpttoward = (TH1F*)fHistos->FindObject("histo23_toward");
848 TH1F * hpttransverse = (TH1F*)fHistos->FindObject("histo23_transverse");
849 TH1F * hptaway = (TH1F*)fHistos->FindObject("histo23_away");
850
851 if ( (hprof1toward) && (hprof1transverse) && (hprof1away) && (hprof2toward) && (hprof2transverse) && (hprof2away) )
852 {
853 for( Int_t part = 0 ; part < fNPart ; part++)
854 {
855 Double_t ptpart = fVectParticle[part]->pt ; // pt of particle
856 if ( ( dphipartljet[part] >=0.) && ( dphipartljet[part] < kPI/3. ) )
857 {
858 hprof1toward->Fill( ptleadjet, fNPart );
859 hprof2toward->Fill( ptleadjet, ptsumevent);
860 hpttoward->Fill( ptpart );
861 }
862 else
863 if ( ( dphipartljet[part] >= (kPI/3.)) && ( dphipartljet[part] < (2.*kPI/3.)) )
864 {
865 hprof1transverse->Fill( ptleadjet, fNPart );
866 hprof2transverse->Fill( ptleadjet, ptsumevent);
867 hpttransverse->Fill( ptpart );
868 }
869 else
870 if ( ( dphipartljet[part] >= ( 2.*kPI/3.)) && ( dphipartljet[part] < kPI ) )
871 {
872 hprof1away->Fill( ptleadjet, fNPart );
873 hprof2away->Fill( ptleadjet, ptsumevent);
874 hptaway->Fill( ptpart );
875 }
876 }
877 }
878
879 delete [] dphipartljet;
880 delete [] zpartljet;
881 delete [] idxpartLJ;
882
883}
884
885//______________________________________________________________________________
886void AliCdfJetFinder::Clean()
887{
888 // CLEANING SECTION
889 for ( Int_t i = 0 ; i < fNPart ; i++ ){
890 delete fVectParticle[i];
891 fVectParticle[i] = 0;
892 }
893 delete [] fVectParticle;fVectParticle = 0;
894
895 for ( Int_t i = 0 ; i < fNJets ; i++ ){
896 delete fVectJet[i];
897 fVectJet[i] = 0;
898 }
899 delete [] fVectJet;fVectJet = 0;
900
901 delete [] fPtArray;fPtArray = 0;
902 delete [] fIdxArray;fIdxArray = 0;
903
904 Reset();
905
906}
907