Fixes to track selection for embedding in cluster task (Bastian), Do no select events...
[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),
8ddcf605 67 fNJets(0),
68 fNPart(0),
7c6659fc 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
8ddcf605 894 delete [] dphipartljet;
895 delete [] zpartljet;
896 delete [] idxpartLJ;
7c6659fc 897}
898
899
900//______________________________________________________________________________
901void AliCdfJetFinder::Clean()
902{
903// CLEANING SECTION
95345ec0 904 for ( Int_t i = 0 ; i < fNPart ; i++ ){
905 delete fVectParticle[i];
906 fVectParticle[i] = 0;
907 }
908 delete [] fVectParticle;fVectParticle = 0;
909
910 for ( Int_t i = 0 ; i < fNJets ; i++ ){
911 delete fVectJet[i];
912 fVectJet[i] = 0;
913 }
914 delete [] fVectJet;fVectJet = 0;
915
916 delete [] fPtArray;fPtArray = 0;
917 delete [] fIdxArray;fIdxArray = 0;
7c6659fc 918
919 Reset();
920}
921
922//______________________________________________________________________________
923void AliCdfJetFinder::FinishRun()
924{
925// do i need this?
926}
927