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