1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
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 //---------------------------------------------------------------------
35 #include <Riostream.h>
39 #include <TClonesArray.h>
40 #include <TLorentzVector.h>
45 #include "AliCdfJetFinder.h"
46 #include "AliCdfJetHeader.h"
47 #include "AliJetReader.h"
48 #include "AliJetReaderHeader.h"
50 #include "AliAODJet.h"
51 #include "AliAODEvent.h"
55 ClassImp ( AliCdfJetFinder )
57 //______________________________________________________________________________
58 AliCdfJetFinder::AliCdfJetFinder():
77 //______________________________________________________________________________
78 AliCdfJetFinder::~AliCdfJetFinder()
84 //______________________________________________________________________________
85 void AliCdfJetFinder::CreateOutputObjects(TList *histos)
87 // Create the list of histograms. Only the list is owned.
90 // gStyle->SetOptStat(11111111);
92 TH1F *h1 = new TH1F ("histo1", "Pt distribution of jets", 200, 0,200); // 1GeV/bin
94 h1->GetXaxis()->SetTitle("P_{T} of jets");
95 h1->GetYaxis()->SetTitle("Number of jets");
96 h1->GetXaxis()->SetTitleColor(1);
97 h1->SetMarkerStyle(kFullCircle);
100 TH1F *h2 = new TH1F ("histo2", "Eta distribution of jets", 240, -1.2,1.2); // 1 unit of rapidity / 100 bin
102 h2->GetXaxis()->SetTitle("Eta of jets");
103 h2->GetYaxis()->SetTitle("Number of jets");
104 h2->GetXaxis()->SetTitleColor(1);
105 h2->SetMarkerStyle(kFullCircle);
108 TH1F *h3 = new TH1F ("histo3", "Phi distribution of jets", 400, -4,4);
110 h3->GetXaxis()->SetTitle("Phi of jets");
111 h3->GetYaxis()->SetTitle("Number of jets");
112 h3->GetXaxis()->SetTitleColor(1);
113 h3->SetMarkerStyle(kFullCircle);
116 TH1F *h4 = new TH1F ("histo4", "Multiplicity of jets", 40, 0,40); // 1 unit of multiplicity /bin
118 h4->GetXaxis()->SetTitle("Particles in jets");
119 h4->GetYaxis()->SetTitle("Number of jets");
120 h4->GetXaxis()->SetTitleColor(1);
121 h4->SetMarkerStyle(kFullCircle);
124 TH1F *h5 = new TH1F ("histo5", "Distribution of jets in events", 100, 0,100);
126 h5->GetXaxis()->SetTitle("Number of jets");
127 h5->GetYaxis()->SetTitle("Number of events");
128 h5->GetXaxis()->SetTitleColor(1);
129 h5->SetMarkerStyle(kFullCircle);
132 TH1F *h6 = new TH1F ("histo6", "Jet1 Charged Multiplicity Distribution", 30, 0,30);
134 h6->GetXaxis()->SetTitle("N_{chg}");
135 h6->GetYaxis()->SetTitle("Number of jets");
136 h6->GetXaxis()->SetTitleColor(1);
137 h6->SetMarkerStyle(kFullCircle);
140 TProfile * h7 = new TProfile ("histo7","N_{chg}(jet1) vs P_{T}(charged jet1)", 200, 0. ,200. , 0.,200. ) ;
142 h7->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
143 h7->GetYaxis()->SetTitle("<N_{chg}(jet1)> in 1 GeV/c bin");
144 h7->GetXaxis()->SetTitleColor(1);
145 h7->SetMarkerStyle(kFullCircle);
148 TH1F *h8 = new TH1F ("histo8", "Charge momentum distribution for leading jet", 120, 0 , 1.2);
150 h8->GetXaxis()->SetTitle("Jets");
151 h8->GetYaxis()->SetTitle("Particle distribution");
152 h8->GetXaxis()->SetTitleColor(1);
153 h8->SetMarkerStyle(kFullCircle);
156 TProfile *h9 = new TProfile ("histo9", "N_{chg} vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 20 );
158 h9->GetXaxis()->SetTitle("#Delta#phi (degrees)");
159 h9->GetYaxis()->SetTitle("<N_{chg}> in 3.6 degree bin");
160 h9->GetXaxis()->SetTitleColor(1);
161 h9->SetMarkerStyle(kFullCircle);
164 TProfile *h10 = new TProfile ("histo10", "P_{T} sum vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 100 );
165 h10->SetStats(kTRUE);
166 h10->GetXaxis()->SetTitle("#Delta#phi (degrees)");
167 h10->GetYaxis()->SetTitle("<P_{T} sum> in 3.6 degree bin");
168 h10->GetXaxis()->SetTitleColor(1);
169 h10->SetMarkerStyle(kFullCircle);
172 TH1F *h11 = new TH1F ("histo11", " \"Transverse\" Pt Distribution ", 70, 0 , 14);
173 h11->SetStats(kTRUE);
174 h11->GetXaxis()->SetTitle("P_{T} (GeV/c)");
175 h11->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
176 h11->GetXaxis()->SetTitleColor(1);
177 h11->SetMarkerStyle(kFullCircle);
180 TH1F *h20 = new TH1F ("histo20", "Distribution of R in leading jet", 400, 0.,4.);
181 h20->SetStats(kTRUE);
182 h20->GetXaxis()->SetTitle("R [formula]");
183 h20->GetYaxis()->SetTitle("dN/dR");
184 h20->GetXaxis()->SetTitleColor(1);
185 h20->SetMarkerStyle(kFullCircle);
188 TProfile * h21 = new TProfile ("histo21","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 30. ) ;
189 h21->SetStats(kTRUE);
190 h21->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
191 h21->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
192 h21->GetXaxis()->SetTitleColor(1);
193 h21->SetMarkerStyle(kFullCircle);
196 TProfile * h22 = new TProfile ("histo22","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
197 h22->SetStats(kTRUE);
198 h22->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
199 h22->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
200 h22->GetXaxis()->SetTitleColor(1);
201 h22->SetMarkerStyle(kFullCircle);
204 TProfile * h21_toward = new TProfile ("histo21_toward","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
205 h21_toward->SetStats(kTRUE);
206 h21_toward->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
207 h21_toward->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
208 h21_toward->GetXaxis()->SetTitleColor(1);
209 h21_toward->SetMarkerStyle(kFullCircle);
210 fHistos->Add(h21_toward);
212 TProfile * h21_transverse = new TProfile ("histo21_transverse","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
213 h21_transverse->SetStats(kTRUE);
214 h21_transverse->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
215 h21_transverse->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
216 h21_transverse->GetXaxis()->SetTitleColor(1);
217 h21_transverse->SetMarkerStyle(kFullCircle);
218 fHistos->Add(h21_transverse);
220 TProfile * h21_away = new TProfile ("histo21_away","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
221 h21_away->SetStats(kTRUE);
222 h21_away->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
223 h21_away->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
224 h21_away->GetXaxis()->SetTitleColor(1);
225 h21_away->SetMarkerStyle(kFullCircle);
226 fHistos->Add(h21_away);
228 TProfile * h22_toward = new TProfile ("histo22_toward","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
229 h22_toward->SetStats(kTRUE);
230 h22_toward->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
231 h22_toward->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
232 h22_toward->GetXaxis()->SetTitleColor(1);
233 h22_toward->SetMarkerStyle(kFullCircle);
234 fHistos->Add(h22_toward);
236 TProfile * h22_transverse = new TProfile ("histo22_transverse","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
237 h22_transverse->SetStats(kTRUE);
238 h22_transverse->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
239 h22_transverse->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
240 h22_transverse->GetXaxis()->SetTitleColor(1);
241 h22_transverse->SetMarkerStyle(kFullCircle);
242 fHistos->Add(h22_transverse);
244 TProfile * h22_away = new TProfile ("histo22_away","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
245 h22_away->SetStats(kTRUE);
246 h22_away->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
247 h22_away->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
248 h22_away->GetXaxis()->SetTitleColor(1);
249 h22_away->SetMarkerStyle(kFullCircle);
250 fHistos->Add(h22_away);
252 TH1F *h23_toward = new TH1F ("histo23_toward","'Toward' Pt Distribution of charged particles", 200, 0., 14.);
253 h23_toward->SetStats(kTRUE);
254 h23_toward->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
255 h23_toward->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
256 h23_toward->GetXaxis()->SetTitleColor(1);
257 h23_toward->SetMarkerStyle(kFullCircle);
258 fHistos->Add(h23_toward);
260 TH1F *h23_transverse = new TH1F ("histo23_transverse","'Transverse' Pt Distribution of charged particles", 200, 0., 14.);
261 h23_transverse->SetStats(kTRUE);
262 h23_transverse->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
263 h23_transverse->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
264 h23_transverse->GetXaxis()->SetTitleColor(1);
265 h23_transverse->SetMarkerStyle(kFullCircle);
266 fHistos->Add(h23_transverse);
268 TH1F *h23_away = new TH1F ("histo23_away","'Away' Pt Distribution of charged particles", 200, 0., 14.);
269 h23_away->SetStats(kTRUE);
270 h23_away->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
271 h23_away->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
272 h23_away->GetXaxis()->SetTitleColor(1);
273 h23_away->SetMarkerStyle(kFullCircle);
274 fHistos->Add(h23_away);
276 TProfile * h24 = new TProfile ("histo24","Jet1 Size vs P_{T}(charged jet1)", 200, 0., 50. , 0., 0.5) ;
277 h24->SetStats(kTRUE);
278 h24->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
279 h24->GetYaxis()->SetTitle("<R(chgjet1)> in 1 GeV/c bin");
280 h24->GetXaxis()->SetTitleColor(1);
281 h24->SetMarkerStyle(kFullCircle);
284 TProfile * h25 = new TProfile ("histo25","Jet1 Size vs P_{T}(charged jet1)", 200, 0., 50. , 0., 0.5) ;
285 h25->SetStats(kTRUE);
286 h25->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
287 h25->GetYaxis()->SetTitle("<R(chgjet1)> in 1 GeV/c bin");
288 h25->GetXaxis()->SetTitleColor(1);
289 h25->SetMarkerStyle(kFullCircle);
292 TProfile *h26 = new TProfile ("histo26", "N_{chg} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0., 0.8);
293 h26->SetStats(kTRUE);
294 h26->GetXaxis()->SetTitle("Distance R");
295 h26->GetYaxis()->SetTitle("<N_{chg}> in 0.02 bin");
296 h26->GetXaxis()->SetTitleColor(1);
297 h26->SetMarkerStyle(kFullCircle);
300 TProfile *h27 = new TProfile ("histo27", "N_{chg} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0., 0.8);
301 h27->SetStats(kTRUE);
302 h27->GetXaxis()->SetTitle("Distance R");
303 h27->GetYaxis()->SetTitle("<N_{chg}> in 0.02 bin");
304 h27->GetXaxis()->SetTitleColor(1);
305 h27->SetMarkerStyle(kFullCircle);
308 TProfile *h28 = new TProfile ("histo28", "PT_{sum} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0.01, 10.);
309 h28->SetStats(kTRUE);
310 h28->GetXaxis()->SetTitle("Distance R");
311 h28->GetYaxis()->SetTitle("<PT_{sum} (GeV/c)> in 0.02 bin");
312 h28->GetXaxis()->SetTitleColor(1);
313 h28->SetMarkerStyle(kFullCircle);
316 TProfile *h29 = new TProfile ("histo29", "PT_{sum} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0.01, 10.);
317 h29->SetStats(kTRUE);
318 h29->GetXaxis()->SetTitle("Distance R");
319 h29->GetYaxis()->SetTitle("<PT_{sum} (GeV/c)> in 0.02 bin");
320 h29->GetXaxis()->SetTitleColor(1);
321 h29->SetMarkerStyle(kFullCircle);
327 //______________________________________________________________________________
330 * Order all charged particles according to their PT.
331 * Start with the highest PT particle and include in the "jet" all particles within the "radius" R = 0.7
332 (considering each particle in the order of decreasing PT and recalculating the centroid of the jet after
333 each new particle is added to the jet).
334 * Go to the next highest PT particle (not already included in a jet) and include in the "jet" all particles
335 (not already included in a jet) within the radius R =0.7.
336 * Continue until all particles are in a "jet".
339 //______________________________________________________________________________
341 //______________________________________________________________________________
342 void AliCdfJetFinder::FindJets()
344 AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;
348 fDebug = header->IsDebugCDF();
349 fAODwrite = header->IsAODwrite() ; // write jets to AOD
350 fAODtracksWrite = header->IsAODtracksWrite() ; // write jet tracks to AOD
351 fRadius = header->GetRadius(); // get Radius from jet finder header
352 fMinJetParticles = header->GetMinPartJet (); // get minimum multiplicity of an jet
353 fJetPtCut = header->GetJetPtCut (); // get minimum of jet pt
356 { cout << "Header not found" << endl; return; }
358 // temporary until the other problems are resolved
363 fFromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
364 if (fFromAod) { fRefArr = fReader->GetReferences(); }
367 fFromAod = 0 ; // disable for the moment ; only ESD reading for now
371 if (!fNPart) { if (fDebug) {cout << "entries = 0 ; Event empty !!!" << endl ;} return; } // if event empty then exit
375 ComputeConesWeight();
377 if (fAODwrite) { cout << "Writing AOD" << endl ; WriteJets(); }
385 //______________________________________________________________________________
386 void AliCdfJetFinder::InitData()
389 TClonesArray * vectArray = fReader->GetMomentumArray() ;
390 if ( vectArray == 0 ) { cout << "Could not get the momentum array" << endl; return; }
392 fNPart = vectArray->GetEntries() ; // n particles in this event;
394 if ( fNPart == 0 ) { return ; } // if event empty then exit
396 fJets->SetNinput ( fNPart ) ; // number of input objects
398 fVectParticle = new varContainer* [fNPart]; // container for Particles
399 fVectJet = new varContainer* [fNPart]; // container for Jets
401 fPtArray = new Double_t [fNPart] ; // momentum array
402 fIdxArray = new Int_t [fNPart] ; // index array of sorted pts
404 // initialisation of momentum and index arrays
405 for ( Int_t i = 0 ; i < fNPart ; i++ )
406 {// SORTING STEP :: fPtArray with data from TClonesArray of TLorentzVector
407 TLorentzVector * lv = (TLorentzVector*) vectArray->At(i);
409 // INITIALISATION of local arrays for temporary storage
410 varContainer *aParticle = new varContainer;
411 aParticle->pt = lv->Pt();
412 aParticle->eta = lv->Eta();
413 aParticle->phi = DeltaPhiNorm ( lv->Phi() ); // normalize to -pi,pi
414 aParticle->njet = -999;
416 fVectParticle[i] = aParticle; // vector of Particles
418 // initializing arrays
419 fIdxArray [i] = -999 ;
420 fPtArray [i] = aParticle->pt ;
423 TMath::Sort ( fNPart, fPtArray, fIdxArray ) ; // get a sorted array of indexes with TClonesArray.Size()
428 //______________________________________________________________________________
429 void AliCdfJetFinder::FindCones()
431 Double_t pt_seed = 0. , eta_seed = 0. , phi_seed = 0. ; // leading particle params
432 Double_t pt_tmp = 0. , eta_tmp = 0. , phi_tmp = 0. ; // temporary variables to be used in various calculations
433 Double_t deta = 0. , dphi = 0. , d_computed = 0. ;
436 fNJets = -1 ; // n jets in this event
437 Int_t idxPtSort = -1 ; // index of array of sorted pt indexes
439 if (fDebug) { cout << "\n\n\n\n\n\n------------------\nBegin Event Analysis\n------------------\n\n" << endl ;}
441 cout << "fNPart = " << fNPart << endl;
443 TBits lkup_table ( fNPart ) ; // bit container ; 1-to-1 corespondence with fIdxArray
445 while ( lkup_table.CountBits() != (UInt_t)fNPart )
446 { // loop over particles in event until all flags are set
447 UInt_t first_non_flagged = lkup_table.FirstNullBit() ; // set the index to the first NON flagged bit ; less conditions
449 cout << "\n\nfirst_non_flagged : " << first_non_flagged << endl;
451 ++fNJets; // incrementing the jet counter
452 if (fDebug) { printf("JET %d \n", fNJets); }
454 pt_seed = 0. ; eta_seed = 0. ; phi_seed = 0. ; // reseting leading particle params
456 for ( UInt_t i_part = first_non_flagged ; i_part < (UInt_t)fNPart ; i_part++ )
457 {// iteration over particles in event
458 // the loop is done over sorted array of pt
459 idxPtSort = fIdxArray[i_part] ; // index of particle ! fIdxArray is an index list pt sorted
461 if ( lkup_table.TestBitNumber(i_part) ) { continue; } // if 4vector is already flagged skip it
463 //init computed and used vars
464 pt_tmp = 0. ; eta_tmp = 0. ; phi_tmp = 0. ;
465 deta = 0. ; dphi = 0. ; d_computed = 0. ; in_jet = 0 ;
467 //taking info from fVectParticle ;
468 pt_tmp = fVectParticle[idxPtSort]->pt ;
469 eta_tmp = fVectParticle[idxPtSort]->eta ;
470 phi_tmp = fVectParticle[idxPtSort]->phi ;
472 if ( i_part == first_non_flagged )
473 {// this is first particle in event; leading particle
474 // begin the search around this particle in a fRadius
477 pt_seed = pt_tmp ; eta_seed = eta_tmp ; phi_seed = phi_tmp ; // seeding the jet with first particle idxPtSort
479 lkup_table.SetBitNumber ( i_part ) ; // flag the index of particle in lkup_table
480 fVectParticle[idxPtSort]->njet = fNJets ; // associate particle with current jet number
482 if (fDebug) { printf("\nLeading particle :: particle index = %d ; at sorted index %d ; in jet %d \n", idxPtSort, i_part, fNJets); }
483 if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pt_tmp, eta_tmp, phi_tmp) ; }
484 if (fDebug) { lkup_table.Print() ;}
486 continue ; // skip to next particle
489 // condition to be in jet
490 deta = eta_tmp - eta_seed ;
491 dphi = DeltaPhiNorm ( phi_tmp - phi_seed ) ; // computing dphi and normalizing to (0,2pi) interval in one step
493 d_computed = Distance (deta, dphi) ; // Distance(fRadius) to (eta,phi) seed
495 in_jet = ( ( fRadius - d_computed ) >= 0.000000001 ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0
498 { // calculus of jet variables
499 lkup_table.SetBitNumber ( i_part ) ; // flag the index of particle in lkup_table
500 fVectParticle[idxPtSort]->njet = fNJets ; // setting in particle list the associated jet
502 if (fDebug) { printf("\njet particle :: particle index = %d ; at sorted index %d ; in jet %d ; found at radius %g ; \n", idxPtSort, i_part, fNJets, d_computed); }
503 if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pt_tmp, eta_tmp, phi_tmp) ; }
504 if (fDebug) { lkup_table.Print() ;}
506 continue ; // skip to next particle
510 // end of iteration over event; one jet definition of content ; jet parameters to be computed later
515 //______________________________________________________________________________
516 void AliCdfJetFinder::ComputeConesWeight()
519 CALCULUS OF JETS Pt, Eta and Phi (centre of weight in (eta,phi) plane)
521 // rescan the vector of particles by identify them by asociate jet number for computing of weight centre
522 // we know : fNJets = the number of jets
524 Double_t pt_jet , eta_jet , phi_jet ; Int_t npartJet ;
525 Double_t pt_tmp = 0. , eta_tmp = 0. , phi_tmp = 0. ; // temporary variables to be used in various calculations
526 Int_t idxPtSort = -999 ; // index of array of sorted pt indexes
528 for( Int_t jet = 0 ; jet < fNJets ; jet++ )
530 if (fDebug) { printf("\n\n--- Computing weight of Jet %d \n", jet ); }
531 npartJet = 0 ; pt_jet = 0. ; eta_jet = 0. ; phi_jet = 0. ; // reset variables for a new computation
533 for ( Int_t i_part = 0 ; i_part < fNPart ; i_part++ )
534 {// iteration over particles in event
535 // the loop is done over sorted array of pt
536 idxPtSort = fIdxArray[i_part] ; // index of particle ! fIdxArray is an index list pt sorted
538 if ( fVectParticle[idxPtSort]->njet == jet )
540 ++npartJet; // incrementing the counter of jet particles
542 //taking info from fVectParticle ;
543 pt_tmp = fVectParticle[idxPtSort]->pt ;
544 eta_tmp = fVectParticle[idxPtSort]->eta ;
545 phi_tmp = DeltaPhiNorm ( fVectParticle[idxPtSort]->phi ) ;
548 eta_jet = ( (pt_jet*eta_jet) + (pt_tmp*eta_tmp) )/(pt_jet + pt_tmp) ;
549 phi_jet = ( (pt_jet*phi_jet) + (pt_tmp*phi_tmp) )/(pt_jet + pt_tmp) ;
552 // add a particle and recalculation of centroid
554 // end of 1 jet computation
556 varContainer *aJet = new varContainer; // Jet container
557 aJet->pt = pt_jet; aJet->eta = eta_jet; aJet->phi = phi_jet; aJet->njet = npartJet; // setting jet vars in container
558 fVectJet[jet] = aJet; // store the number of the jet(fNJets) and increment afterwards
560 if (fDebug) { printf ("=== current jet %d : npartjet= %d ; pt_jet= %g ; eta_jet = %g ; phi_jet = %g \n\n\n",
561 jet, npartJet, pt_jet, eta_jet, phi_jet ) ; }
569 //______________________________________________________________________________
570 void AliCdfJetFinder::WriteJets()
571 { // Writing AOD jets and AOD tracks
573 /* for( Int_t jet_nr = 0 ; jet_nr < fNJets ; jet_nr++ )
575 Double_t pt = 0., eta = 0., phi = 0., // jet variables
576 px = 0., py = 0., pz = 0., en = 0.; // convert to 4-vector
577 pt = fVectJet[ jet_nr ]->pt ; // pt of jet
578 eta = fVectJet[ jet_nr ]->eta ; // eta of jet
579 phi = fVectJet[ jet_nr ]->phi ; // phi of jet
581 px = pt * TMath::Cos ( phi ) ;
582 py = pt * TMath::Sin ( phi ) ;
583 pz = pt / TMath::Tan ( 2.0 * TMath::ATan ( TMath::Exp ( -eta ) ) ) ;
584 en = TMath::Sqrt ( px * px + py * py + pz * pz );
586 AliAODJet jet (px, py, pz, en);
589 if (fDebug) jet.Print("");
591 if (fromAod && AODtracksWrite)
593 for ( Int_t jet_track = 0; jet_track < fNPart; jet_track++ )
594 { if ( fVectParticle[jet_track]->njet == jet_nr ) { jet.AddTrack(refs->At(jet_track)) ; } }
596 // tracks REFs written in AOD
600 //jets vector parsed and written to AOD
605 //______________________________________________________________________________
606 void AliCdfJetFinder::AnalizeJets()
609 //persistent pointer to histo20
610 TH1F *h_r = (TH1F*)fHistos->FindObject("histo20");
612 Int_t *jets_pt_idx = 0; // sorted array of jets pt
613 Double_t *jets_pt = 0; // array of jets pts
614 Int_t leading_jet_index = -1 ; // index of leading jet from fVectJet
615 Int_t part_leadjet = 0 ; // number of particles in leading jet
616 Double_t pt_leadjet = 0. ; // pt of leading jet
617 Double_t eta_leadjet = 0. ; // eta of leading jet
618 Double_t phi_leadjet = 0. ; // phi of leading jet
620 jets_pt_idx = new Int_t [fNJets] ;
621 jets_pt = new Double_t [fNJets] ;
623 //________________________________________________________________________________________
624 // Jet sorting and finding the leading jet that coresponds to cuts in pt and multiplicity
625 //________________________________________________________________________________________
627 // filing the idx_ptjets array
628 if (fDebug) printf("List of unsorted jets:\n");
629 for( Int_t i = 0 ; i < fNJets ; i++ )
631 jets_pt_idx [i] = 0 ;
632 jets_pt [i] = fVectJet[i]->pt ;
633 if (fDebug) { cout << " jet found: " << i << " npartjet=" << fVectJet[i]->njet << " ; jets_pt = " << jets_pt[i] << endl; }
636 TMath::Sort ( fNJets, jets_pt , jets_pt_idx ) ; // sorting pt of jets
638 // selection of leading jet
639 // looping over jets searching for __first__ one that coresponds to cuts
640 for( Int_t i = 0 ; i < fNJets ; i++ )
642 if ( ( fVectJet[ jets_pt_idx[i] ]->njet >= fMinJetParticles ) && ( fVectJet[ jets_pt_idx[i] ]->pt >= fJetPtCut ) )
644 leading_jet_index = jets_pt_idx[i] ;
645 part_leadjet = fVectJet[ leading_jet_index ]->njet ; // number of particles in leading jet
646 pt_leadjet = fVectJet[ leading_jet_index ]->pt ; // pt of leading jet
647 eta_leadjet = fVectJet[ leading_jet_index ]->eta ; // eta of leading jet
648 phi_leadjet = fVectJet[ leading_jet_index ]->phi ; // phi of leading jet
651 { printf("Leading jet %d : npart= %d ; pt= %g ; eta = %g ; phi = %g \n", leading_jet_index, part_leadjet, pt_leadjet, eta_leadjet, phi_leadjet ); }
656 // end of selection of leading jet
660 //////////////////////////////////////////////////
661 //// Computing of values used in histograms
662 //////////////////////////////////////////////////
664 //___________________________________________________________________________
665 // pt_sum of all particles in event
666 //___________________________________________________________________________
667 cout << "Computing sum of pt in event" << endl ;
668 Double_t pt_sum_event = 0.;
669 for ( Int_t i = 0 ; i< fNPart ; i++ ) { pt_sum_event += fVectParticle[i]->pt ; }
670 printf ("Sum of all Pt in event : pt_sum_event = %g", pt_sum_event) ;
672 //___________________________________________________________________________
673 // Filling an array with indexes of leading jet particles
674 //___________________________________________________________________________
675 Int_t * idx_leadjet_part = new Int_t [part_leadjet] ;
676 Int_t counter_part_lead_jet = 0;
678 cout << "Filling an array with indexes of leading jet particles" << endl;
680 for( Int_t i = 0 ; i < fNPart ; i++ )
682 if ( fVectParticle[i]->njet == leading_jet_index )
683 { idx_leadjet_part[counter_part_lead_jet++] = i ; }
686 if ( (counter_part_lead_jet-1) > part_leadjet ) { cout << " Counter_part_lead_jet > part_leadjet !!!!" << endl;}
689 //___________________________________________________________________________
690 // Calculus of part distribution in leading jet
691 //___________________________________________________________________________
693 Double_t *z_part_ljet = new Double_t [ part_leadjet ] ; // array of z of particles in leading jet
695 cout << "Entering loop of calculus of part distribution in leading jet" << endl ;
697 for( Int_t j = 0 ; j < part_leadjet ; j++ )
699 Double_t z_j = fVectParticle[idx_leadjet_part[j]]->pt ;
700 z = z_j / pt_leadjet ;
701 z_part_ljet [j] = z ;
702 cout << "idx_leadjet_part[j] = " << idx_leadjet_part[j]
703 << " p of particle = " << z_j
704 << " pt lead jet = " << pt_leadjet
705 << " Z = " << z << endl;
709 //___________________________________________________________________________
710 // array of delta phi's between phi of particles and leading jet phi
711 //___________________________________________________________________________
712 cout << "array of delta phi's between phi of particles and leading jet phi" << endl;
713 Double_t dphi_partLJ = 0. ;
714 Double_t *dphi_part_ljet = new Double_t [fNPart];
715 for( Int_t part = 0 ; part < fNPart ; part++ )
717 dphi_partLJ = fVectParticle[part]->phi - phi_leadjet ;
718 dphi_partLJ = DeltaPhiNorm (dphi_partLJ) ; // restrict the delta phi to (0,pi) interval
719 dphi_part_ljet [part] = dphi_partLJ ;
720 printf("part= %d ; dphi_partLJ = %g \n", part, dphi_partLJ );
724 //______________________________________________________________________________
725 // Pt distribution for all particles
726 //______________________________________________________________________________
727 TH1F * h_pt = (TH1F*)fHistos->FindObject("histo11");
728 if ( h_pt ) { for ( Int_t i = 0 ; i < fNPart ; i++ ) { h_pt->Fill( fVectParticle[i]->pt ); } }
730 //___________________________________________________________________________
731 // Recomputing of radius of particles in leading jet
732 //___________________________________________________________________________
733 if (fDebug) { printf(" Searching particles with jet index %d\n", leading_jet_index); }
735 Double_t ddeta = 0. , ddphi = 0. , r_part = 0. ;
737 for( Int_t j = 0 ; j < part_leadjet ; j++ )
739 ddeta = eta_leadjet - fVectParticle[idx_leadjet_part[j]]->eta;
741 Double_t phi_tmp = fVectParticle[idx_leadjet_part[j]]->phi ;
742 phi_tmp = DeltaPhiNorm (phi_tmp);
744 ddphi = DeltaPhiNorm ( phi_leadjet - phi_tmp ) ; // restrict the delta phi to (-pi,pi) interval
746 r_part = Distance (ddeta, ddphi) ;
748 printf ("Particle %d with Re-Computed radius = %f ", idx_leadjet_part[j], r_part) ;
749 if ( (r_part - fRadius) >= 0.00000001 )
750 { printf (" bigger than selected radius of %f\n", fRadius ); }
754 if (h_r) h_r->Fill(r_part);
760 //_______________________________________________________________________
761 // Computing of radius that contain 80% of Leading Jet ( PT and multiplicity )
762 //_______________________________________________________________________
763 Double_t core_part_leadjet = 0.8 * part_leadjet ;
764 Double_t core_pt_leadjet = 0.8 * pt_leadjet ;
765 Int_t counter_core_part = 0 ;
766 Double_t counter_core_pt = 0. ;
767 Int_t sorted_index = -1 ;
769 TProfile * h_prof_24 = (TProfile*)fHistos->FindObject("histo24");
770 TProfile * h_prof_25 = (TProfile*)fHistos->FindObject("histo25");
772 TProfile * h_prof_26 = (TProfile*)fHistos->FindObject("histo26");
773 TProfile * h_prof_27 = (TProfile*)fHistos->FindObject("histo27");
774 TProfile * h_prof_28 = (TProfile*)fHistos->FindObject("histo28");
775 TProfile * h_prof_29 = (TProfile*)fHistos->FindObject("histo29");
778 if ((h_prof_24) && (h_prof_25) && (h_prof_26) && (h_prof_27) && (h_prof_28) && (h_prof_29) )
780 for( Int_t part = 0 ; part < fNPart ; part++ )
782 Double_t pt_tmp = 0. ; Double_t eta_tmp = 0. ; Double_t phi_tmp = 0. ; // temporary variables
783 Double_t d_part = 0. ;
784 sorted_index = fIdxArray[part] ;
786 if ( fVectParticle [ sorted_index ]->njet == leading_jet_index )
788 pt_tmp = fVectParticle[sorted_index]->pt ;
789 eta_tmp = fVectParticle[sorted_index]->eta ;
790 phi_tmp = fVectParticle[sorted_index]->phi ;
792 ++counter_core_part ;
793 counter_core_pt += pt_tmp ;
795 d_part = Distance ( eta_leadjet - eta_tmp, phi_leadjet - phi_tmp ) ;
797 if ( counter_core_part <= core_part_leadjet ) { h_prof_24->Fill(pt_leadjet, d_part); }
798 if ( counter_core_pt <= core_pt_leadjet ) { h_prof_25->Fill(pt_leadjet, d_part); }
800 if (pt_leadjet > 5.) { h_prof_26->Fill(d_part, counter_core_part); h_prof_28->Fill(d_part, counter_core_pt); }
801 if (pt_leadjet > 30.) { h_prof_27->Fill(d_part, counter_core_part); h_prof_29->Fill(d_part, counter_core_pt); }
816 TH1F *h_jet_pt = (TH1F*)fHistos->FindObject("histo1");
817 TH1F *h_jet_eta = (TH1F*)fHistos->FindObject("histo2");
818 TH1F *h_jet_phi = (TH1F*)fHistos->FindObject("histo3");
819 TH1F *h_jet_njet = (TH1F*)fHistos->FindObject("histo4");
821 for( Int_t jet = 0 ; jet < fNJets ; jet++ )
823 if (h_jet_pt) h_jet_pt ->Fill ( fVectJet[ jet ]->pt ) ;
824 if (h_jet_eta) h_jet_eta ->Fill ( fVectJet[ jet ]->eta ) ;
825 if (h_jet_phi) h_jet_phi ->Fill ( fVectJet[ jet ]->phi ) ;
826 if (h_jet_njet) h_jet_njet ->Fill ( fVectJet[ jet ]->njet ) ;
829 TH1F *h_jets = (TH1F*)fHistos->FindObject("histo5");
830 if (h_jets) h_jets->Fill(fNJets);
832 TH1F *h_leadpart = (TH1F*)fHistos->FindObject("histo6");
833 if (h_leadpart) h_leadpart->Fill(part_leadjet);
835 TProfile * h_prof = (TProfile*)fHistos->FindObject("histo7");
836 if (h_prof) h_prof->Fill(pt_leadjet,part_leadjet);
838 TH1F *h_MD = (TH1F*)fHistos->FindObject("histo8");
839 for( Int_t k = 0 ; k < part_leadjet ; k++)
840 { h_MD->Fill( z_part_ljet[k] ); }
842 TProfile * h_phi = (TProfile*)fHistos->FindObject("histo9");
843 for( Int_t k = 0 ; k < part_leadjet ; k++)
844 { h_phi->Fill( TMath::RadToDeg() * dphi_part_ljet [k] , fNPart ) ; }
846 TProfile * h_tpd = (TProfile*)fHistos->FindObject("histo10");
847 for( Int_t k = 0 ; k < fNPart ; k++)
848 { h_tpd->Fill( TMath::RadToDeg() * dphi_part_ljet [k] , pt_sum_event ) ; }
851 TProfile * h_prof1 = (TProfile*)fHistos->FindObject("histo21");
852 if (h_prof1) h_prof1->Fill(pt_leadjet, fNPart);
854 TProfile * h_prof2 = (TProfile*)fHistos->FindObject("histo22");
855 if (h_prof2) h_prof2->Fill(pt_leadjet, pt_sum_event);
857 TProfile * h_prof1_toward = (TProfile*)fHistos->FindObject("histo21_toward");
858 TProfile * h_prof1_transverse = (TProfile*)fHistos->FindObject("histo21_transverse");
859 TProfile * h_prof1_away = (TProfile*)fHistos->FindObject("histo21_away");
860 TProfile * h_prof2_toward = (TProfile*)fHistos->FindObject("histo22_toward");
861 TProfile * h_prof2_transverse = (TProfile*)fHistos->FindObject("histo22_transverse");
862 TProfile * h_prof2_away = (TProfile*)fHistos->FindObject("histo22_away");
863 TH1F * h_pt_toward = (TH1F*)fHistos->FindObject("histo23_toward");
864 TH1F * h_pt_transverse = (TH1F*)fHistos->FindObject("histo23_transverse");
865 TH1F * h_pt_away = (TH1F*)fHistos->FindObject("histo23_away");
869 if ( (h_prof1_toward) && (h_prof1_transverse) && (h_prof1_away) && (h_prof2_toward) && (h_prof2_transverse) && (h_prof2_away) )
871 for( Int_t part = 0 ; part < fNPart ; part++)
873 Double_t pt_part = fVectParticle[part]->pt ; // pt of particle
874 if ( ( dphi_part_ljet[part] >=0.) && ( dphi_part_ljet[part] < pi/3. ) )
876 h_prof1_toward->Fill( pt_leadjet, fNPart );
877 h_prof2_toward->Fill( pt_leadjet, pt_sum_event);
878 h_pt_toward->Fill( pt_part );
881 if ( ( dphi_part_ljet[part] >= (pi/3.)) && ( dphi_part_ljet[part] < (2.*pi/3.)) )
883 h_prof1_transverse->Fill( pt_leadjet, fNPart );
884 h_prof2_transverse->Fill( pt_leadjet, pt_sum_event);
885 h_pt_transverse->Fill( pt_part );
888 if ( ( dphi_part_ljet[part] >= ( 2.*pi/3.)) && ( dphi_part_ljet[part] < pi ) )
890 h_prof1_away->Fill( pt_leadjet, fNPart );
891 h_prof2_away->Fill( pt_leadjet, pt_sum_event);
892 h_pt_away->Fill( pt_part );
906 //______________________________________________________________________________
907 void AliCdfJetFinder::Clean()
912 delete [] fVectParticle;
917 // if (z_part_ljet) delete [] z_part_ljet;
918 // if (idx_leadjet_part) delete [] idx_leadjet_part;
919 // if (jets_pt_idx) delete [] jets_pt_idx ;
920 // if (jets_pt) delete [] jets_pt ;
922 // vectArray->Delete() ; // TClonesArray of lorentz vectors
933 //______________________________________________________________________________
934 void AliCdfJetFinder::FinishRun()