1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
16 //---------------------------------------------------------------------
\r
19 //---------------------------------------------------------------------
\r
22 #include <Riostream.h>
\r
24 #include <TCanvas.h>
\r
26 #include <TClonesArray.h>
\r
27 #include <TLorentzVector.h>
\r
31 #include <TArrayF.h>
\r
32 #include "AliCdfJetFinder.h"
\r
33 #include "AliCdfJetHeader.h"
\r
34 #include "AliJetReader.h"
\r
35 #include "AliJetReaderHeader.h"
\r
37 #include "AliAODJet.h"
\r
38 #include "AliAODEvent.h"
\r
39 #include "TProfile.h"
\r
41 using namespace std ;
\r
43 struct varContainer // container for Particle variables
\r
44 { // variables of container struct
\r
45 Double_t pt; Double_t eta; Double_t phi;
\r
46 Int_t njet; // if jets are stored in varContainer nr_jet is multiplicity of jet
\r
49 ClassImp ( AliCdfJetFinder )
\r
51 //______________________________________________________________________________
\r
52 AliCdfJetFinder::AliCdfJetFinder():
\r
55 { /* Constructor */ }
\r
57 //______________________________________________________________________________
\r
58 AliCdfJetFinder::~AliCdfJetFinder()
\r
62 cout << "Calling the destructor ... " << endl ;
\r
64 cout << "Destructor called!" << endl;
\r
67 //______________________________________________________________________________
\r
68 void AliCdfJetFinder::CreateOutputObjects(TList *histos)
\r
70 // Create the list of histograms. Only the list is owned.
\r
73 TH1F *h1 = new TH1F ("histo1", "no jets/event", 100, 0,20);
\r
74 h1->SetStats(kTRUE);
\r
75 h1->GetXaxis()->SetTitle("N_{jets}");
\r
76 h1->GetYaxis()->SetTitle("#frac{dN}{dN_{jets}}");
\r
77 h1->GetXaxis()->SetTitleColor(1);
\r
78 h1->SetMarkerStyle(kFullCircle);
\r
81 TH1F *h2 = new TH1F ("histo2", "no part/leading jet", 40, 0,40);
\r
82 h2->SetStats(kTRUE);
\r
83 h2->GetXaxis()->SetTitle("N_{leading}");
\r
84 h2->GetYaxis()->SetTitle("#frac{dN}{dN_{leading}}");
\r
85 h2->GetXaxis()->SetTitleColor(1);
\r
86 h2->SetMarkerStyle(kFullCircle);
\r
89 TProfile * h3 = new TProfile ("histo3","ProfileX of (pt,npart) of leading jets", 25, 0. ,25. , 0.,25. ) ;
\r
90 h3->SetStats(kTRUE);
\r
91 h3->GetXaxis()->SetTitle("P_{t}");
\r
92 h3->GetYaxis()->SetTitle("N_{leading}");
\r
93 h3->GetXaxis()->SetTitleColor(1);
\r
94 h3->SetMarkerStyle(kFullCircle);
\r
97 TH1F *h4 = new TH1F ("histo4", "Charge momentum distribution for leading jet", 500, 0 , 5);
\r
100 TProfile *h5 = new TProfile ("histo5", "ProfileX of N_{charge} vs dphi leading jet", 100 , 0. , 200. , 0 , 100 );
\r
103 TH1F *h6 = new TH1F ("histo6", " \"Transverse\" Pt Distribution ", 70, 0 , 14);
\r
107 //cout << "CreateOutputObjetcs done!" << endl ;
\r
110 //______________________________________________________________________________
\r
111 void AliCdfJetFinder::FindJets()
\r
113 //cout << "Begining of FindJets ..." << endl ;
\r
114 //1) Fill 1 array with momentum and initialise another array for indexes
\r
116 // 2.1) Sort momentum array
\r
117 // 2.2) loop over arrays of TLorentzVectors
\r
118 // 2.3) flag as a possible jet
\r
119 //3) fill AliAODJet list
\r
120 Bool_t debug = kFALSE;
\r
122 TRefArray *refs = 0;
\r
123 Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
\r
124 if (fromAod) { refs = fReader->GetReferences(); }
\r
126 AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;
\r
127 Double_t radius = header->GetRadius(); // get Radius from header
\r
128 // cout << "Radius is : " << radius << endl ;
\r
130 TClonesArray * vectArray = fReader->GetMomentumArray() ;
\r
131 if ( vectArray == 0 ) { cout << "Could not get the momentum array" << endl; return; }
\r
133 Int_t nPart = vectArray->GetEntries() ; // n particles in this event;
\r
134 if ( nPart == 0 ) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; }
\r
136 Int_t nJets = 0; // n jets in this event
\r
138 fJets->SetNinput ( nPart ) ; // number of input objects
\r
140 varContainer **vectParticle = new varContainer* [nPart]; // container for Particles
\r
141 varContainer **vectJet = new varContainer* [nPart]; // container for Jets
\r
143 Double_t *ptArray = new Double_t [nPart] ; // momentum array
\r
144 Int_t *idxArray = new Int_t [nPart] ; // index array of sorted pts
\r
146 // initialisation of momentum and index arrays
\r
147 // cout << "Filling idxArray && momArray for sorting " << endl;
\r
148 for ( Int_t i = 0 ; i < nPart ; i++ )
\r
149 {// SORTING STEP :: ptArray with data from TClonesArray of TLorentzVector
\r
150 TLorentzVector * lv = (TLorentzVector*) vectArray->At(i);
\r
151 // INITIALISATION of local arrays for temporary storage
\r
152 varContainer *aParticle = new varContainer;
\r
153 aParticle->pt = lv->Pt();
\r
154 aParticle->eta = lv->Eta();
\r
155 aParticle->phi = lv->Phi();
\r
156 aParticle->njet = -1;
\r
157 vectParticle[i] = aParticle;
\r
159 // initializing arrays
\r
161 ptArray [i] = aParticle->pt ;
\r
164 TMath::Sort ( nPart, ptArray, idxArray ) ; // get a sorted array of indexes with size = size of TClonesArray
\r
166 Double_t pt_jet = 0. ,
\r
170 sum_phi = 0. ; // jet variables
\r
172 Int_t idxPtSort = 0 , // index of array of sorted pt indexes
\r
173 npartJet = 0 ; // number of particles in curent jet
\r
175 TBits lkup_table ( nPart ) ; // bit container of size npart
\r
177 while ( lkup_table.CountBits(0) != (UInt_t)nPart )
\r
178 { // loop over particles in event until all flags are set
\r
179 UInt_t first_non_flagged = lkup_table.FirstNullBit(0) ; // set the index to the first NON flagged bit ; less conditions
\r
181 npartJet = 0 ; // reseting number of particles in jet
\r
188 for ( UInt_t i_part = first_non_flagged ; i_part < (UInt_t)nPart ; i_part++ )
\r
189 {// iteration over particles in event
\r
190 // reseting variables
\r
192 // the loop is done over sorted array of pt
\r
193 idxPtSort = idxArray[i_part] ; // index entry of TLorentzVector from vect_arr
\r
195 if ( lkup_table.TestBitNumber(i_part) == 1 ) { continue; } // if 4vector is already flagged skip it
\r
197 //taking info from vectParticle ;
\r
198 Double_t pt_tmp = vectParticle[idxPtSort]->pt ;
\r
199 Double_t eta_tmp = vectParticle[idxPtSort]->eta ;
\r
200 Double_t phi_tmp = vectParticle[idxPtSort]->phi ;
\r
201 if (debug) printf(" particle %d: pt=%g\n", i_part, pt_tmp);
\r
203 // all angles are expressed in rad
\r
205 if ( TMath::Abs(eta_tmp) > 0.9 )
\r
207 lkup_table.SetBitNumber ( i_part, 1 ) ; // mark particle as used to be skipped
\r
212 if ( i_part == first_non_flagged )
\r
213 {// this is first particle in event; leading particle
\r
214 // initialise the jet variables with leading particle numbers
\r
216 npartJet++; // incrementing counter of particles in jet (nparJet == 1) (leading particle)
\r
218 eta_jet = eta_tmp ;
\r
219 phi_jet = phi_tmp ;
\r
220 if (debug) printf(" first part in jet: npartjet=%d pt_jet=%g\n", npartJet, pt_jet);
\r
222 lkup_table.SetBitNumber ( i_part, 1 ) ; // flag the 4vector
\r
224 vectParticle[idxPtSort]->njet = nJets ; // associate particle with current jet
\r
226 continue ; // skip to next particle
\r
229 // condition to be in jet
\r
230 Double_t deta = eta_jet - eta_tmp ;
\r
231 Double_t dphi = phi_jet - phi_tmp ;
\r
233 if ( dphi < -TMath::Pi() ) { dphi = -dphi - 2.0 * TMath::Pi() ; }
\r
234 if ( dphi > TMath::Pi() ) { dphi = -dphi + 2.0 * TMath::Pi() ; }
\r
236 Double_t r_computed = TMath::Sqrt ( deta * deta + dphi * dphi ) ;
\r
238 Bool_t in_jet = ( r_computed <= radius ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0
\r
241 { // calculus of jet variables
\r
242 npartJet++; // incrementing counter of particles in jet
\r
245 sum_eta += (pt_tmp * eta_tmp) ;
\r
246 sum_phi += (pt_tmp * phi_tmp) ;
\r
247 eta_jet = sum_eta / pt_jet ;
\r
248 phi_jet = sum_phi / pt_jet ;
\r
250 lkup_table.SetBitNumber ( i_part, 1 ) ; // flag the 4vector
\r
251 vectParticle[idxPtSort]->njet = nJets ; // setting in particle list the associated jet
\r
252 if (debug) printf(" particle added to jet: npartjet=%d pt_jet=%g\n", npartJet, pt_jet);
\r
253 continue ; // skip to next particle
\r
257 // end of iteration over event; one jet definition
\r
259 varContainer *aJet = new varContainer;
\r
261 aJet->eta = eta_jet;
\r
262 aJet->phi = phi_jet;
\r
263 aJet->njet = npartJet;
\r
264 vectJet[nJets++] = aJet; // store the jet (from jet 0 to nJets-1 )
\r
265 if (debug) printf("=== current jet: npartjet=%d ptjet=%g\n", npartJet,pt_jet);
\r
267 // writing aod information
\r
268 Double_t px = 0., py = 0., pz = 0., en = 0.; // convert to 4-vector
\r
269 px = pt_jet * TMath::Cos ( phi_jet ) ;
\r
270 py = pt_jet * TMath::Sin ( phi_jet ) ;
\r
271 pz = pt_jet / TMath::Tan ( 2.0 * TMath::ATan ( TMath::Exp ( -eta_jet ) ) ) ;
\r
272 en = TMath::Sqrt ( px * px + py * py + pz * pz );
\r
274 if (npartJet<2) continue; // do not add jets with less than 2 particles
\r
275 if (debug) cout << "Jet 4vect : " << "px = " << px << " ; py = " << py << " ; pz = " << pz << " ; E = " << en << endl;
\r
277 AliAODJet jet (px, py, pz, en);
\r
278 // cout << "Printing jet " << endl;
\r
279 if (debug) jet.Print("");
\r
282 for (Int_t parts = 0; parts < nPart; parts++ ) {
\r
283 if (vectParticle[parts]->njet == nJets-1)
\r
284 jet.AddTrack(refs->At(parts));
\r
288 // cout << "Adding jet ... " ;
\r
290 // cout << "added \n" << endl;
\r
293 // end of while loop over particles ; ends when all particles were flagged as used and all jets defined
\r
295 /////////////////////////////////
\r
296 //// END OF EVENT PARSING ///
\r
297 /////////////////////////////////
\r
299 Int_t *jets_pt_idx = 0; // sorted array of jets pt
\r
300 Double_t *jets_pt = 0; // array of jets pts
\r
301 Int_t leading_jet_index = -1 ; // index of leading jet from vectJet
\r
302 Int_t part_leadjet = 0 ; // number of particles in leading jet
\r
303 Double_t pt_leadjet = 0. ; // pt of leading jet
\r
304 Double_t eta_leadjet = 0. ; // eta of leading jet
\r
305 Double_t phi_leadjet = 0. ; // phi of leading jet
\r
306 Int_t * idx_leadjet_part = 0;
\r
309 { // if there is at least one jet in event
\r
310 jets_pt_idx = new Int_t [nJets] ;
\r
311 jets_pt = new Double_t [nJets] ;
\r
313 // filing the idx_ptjets array
\r
314 if (debug) printf("List of unsorted jets:\n");
\r
315 for( Int_t i = 0 ; i < nJets ; i++ )
\r
317 jets_pt_idx [i] = 0 ;
\r
318 jets_pt [i] = vectJet[i]->pt ;
\r
319 if (debug) cout << " jet found: " << i << " npartjet=" << vectJet[i]->njet << " jets_pt[i]= " << jets_pt [i] << endl;
\r
321 TMath::Sort ( nJets, jets_pt , jets_pt_idx ) ; // sorting pt of jets
\r
323 // selection of leading jet with nr of particles > 1
\r
324 for( Int_t i = 0 ; i < nJets ; i++ )
\r
326 if ( vectJet[ jets_pt_idx[i] ]->njet > 1 )
\r
328 leading_jet_index = jets_pt_idx[i] ;
\r
329 part_leadjet = vectJet[ leading_jet_index ]->njet ; // number of particles in leading jet
\r
330 pt_leadjet = vectJet[ leading_jet_index ]->pt ; // pt of leading jet
\r
331 eta_leadjet = vectJet[ leading_jet_index ]->eta ; // eta of leading jet
\r
332 phi_leadjet = vectJet[ leading_jet_index ]->phi ; // phi of leading jet
\r
334 cout << "Leading jet: npart = " << part_leadjet
\r
335 << " ; pt = " << pt_leadjet
\r
336 << " ; phi = " << phi_leadjet
\r
337 << " ; eta = " << eta_leadjet
\r
342 // end of selection of leading jet
\r
344 // Filling an array with indexes of leading jet particles
\r
345 idx_leadjet_part = new Int_t [part_leadjet] ;
\r
347 if (debug) printf(" Searching particles with jet index %d\n", leading_jet_index);
\r
348 for( Int_t i = 0 ; i <nPart ; i++ )
\r
350 if ( vectParticle[i]->njet == leading_jet_index )
\r
352 idx_leadjet_part[counter++] = i ;
\r
353 if (debug) cout << " " << counter-1 << ": index=" << i << " pt=" << vectParticle[i]->pt << endl;
\r
356 if ( (counter-1) > part_leadjet ) { cout << " Counter > part_leadjet !!!!" << endl;}
\r
359 // end of loop over DETECTED jets
\r
363 // Calculus of part distribution in leading jet
\r
365 Double_t *z_part_ljet = new Double_t [ part_leadjet ] ; // array of z of particles in leading jet
\r
366 for( Int_t j = 0 ; j < part_leadjet ; j++ )
\r
368 Double_t z_j = vectParticle[idx_leadjet_part[j]]->pt ;
\r
369 z = z_j / pt_leadjet ;
\r
370 z_part_ljet [j] = z ;
\r
371 cout << "idx_leadjet_part[j] = " << idx_leadjet_part[j]
\r
372 << " p of particle = " << z_j
\r
373 << " pt lead jet = " << pt_leadjet
\r
374 << " Z = " << z << endl;
\r
379 Double_t dphi_partLJ = 0. ;
\r
380 Double_t dphiArray [nPart];
\r
381 for( Int_t part = 0 ; part < nPart ; part++ )
\r
383 dphi_partLJ = phi_leadjet - phi_partT [part] ;
\r
385 if ( dphi_partLJ < -TMath::Pi() ) { dphi_partLJ = -dphi_partLJ - 2.0 * TMath::Pi() ; }
\r
386 if ( dphi_partLJ > TMath::Pi() ) { dphi_partLJ = -dphi_partLJ + 2.0 * TMath::Pi() ; }
\r
388 dphi_partLJ = TMath::Abs(dphi_partLJ) ;
\r
390 dphiArray [part] = dphi_partLJ ;
\r
393 // Fill nr_jet histogram
\r
396 // printf("FILLING histograms\n");
\r
397 TH1F *h_jets = (TH1F*)fHistos->FindObject("histo1");
\r
398 if (h_jets) h_jets->Fill(nJets);
\r
400 // There may be no leading jet with more that 2 particles !
\r
401 TH1F *h_leadpart = (TH1F*)fHistos->FindObject("histo2");
\r
402 if (part_leadjet && h_leadpart) h_leadpart->Fill(part_leadjet);
\r
404 TProfile * h_prof = (TProfile*)fHistos->FindObject("histo3");
\r
405 if (part_leadjet && h_prof) h_prof->Fill(pt_leadjet,part_leadjet);
\r
407 TH1F *h_MD = (TH1F*)fHistos->FindObject("histo4");
\r
408 for( Int_t k = 0 ; k < part_leadjet ; k++)
\r
409 { h_MD->Fill( z_part_ljet[k] ); }
\r
412 // TProfile *h_phi = (TProfile*)fHistos->At(4);
\r
413 // for( Int_t k = 0 ; k < nPart ; k++)
\r
414 // { h_phi->Fill( TMath::RadToDeg() * dphiArray [k] , nPart ) ; }
\r
417 // TH1F *h_tpd = (TH1F*)fHistos->At(5);
\r
418 // for( Int_t k = 0 ; k < nPart ; k++)
\r
420 // if ( (phi_partT[k] > (TMath::Pi()/3.)) && (phi_partT[k] < (2. * TMath::Pi()/3.)) )
\r
421 // { h_tpd->Fill( pt_partT [k] ) ; }
\r
428 // CLEANING SECTION
\r
429 for (Int_t i=0; i<nPart; i++) delete vectParticle[i];
\r
430 delete [] vectParticle;
\r
431 for (Int_t i=0; i<nJets; i++) delete vectJet[i];
\r
434 delete [] idxArray;
\r
436 if (z_part_ljet) delete [] z_part_ljet;
\r
437 if (idx_leadjet_part) delete [] idx_leadjet_part;
\r
438 if (jets_pt_idx) delete [] jets_pt_idx ;
\r
439 if (jets_pt) delete [] jets_pt ;
\r
441 vectArray->Delete() ; // TClonesArray of lorentz vectors
\r
446 //______________________________________________________________________________
\r
447 void AliCdfJetFinder::FinishRun()
\r