]>
Commit | Line | Data |
---|---|---|
98178771 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to 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 | **************************************************************************/ | |
15 | ||
16 | //--------------------------------------------------------------------- | |
17 | // | |
18 | // | |
19 | //--------------------------------------------------------------------- | |
20 | #include <vector> | |
21 | #include <iostream> | |
22 | #include <Riostream.h> | |
23 | #include <TFile.h> | |
24 | #include <TCanvas.h> | |
25 | #include <TROOT.h> | |
26 | #include <TClonesArray.h> | |
27 | #include <TLorentzVector.h> | |
28 | #include <TH1F.h> | |
29 | #include <TH2F.h> | |
30 | #include <TBits.h> | |
31 | #include <TArrayF.h> | |
32 | #include "AliCdfJetFinder.h" | |
33 | #include "AliCdfJetHeader.h" | |
34 | #include "AliJetReader.h" | |
35 | #include "AliJetReaderHeader.h" | |
36 | #include "AliJet.h" | |
37 | #include "AliAODJet.h" | |
38 | #include "AliAODEvent.h" | |
39 | #include "TProfile.h" | |
40 | ||
41 | using namespace std ; | |
42 | ||
43 | struct varContainer // container for Particle variables | |
44 | { // variables of container struct | |
45 | Double_t pt; Double_t eta; Double_t phi; | |
46 | Int_t njet; // if jets are stored in varContainer nr_jet is multiplicity of jet | |
47 | } ; | |
48 | ||
49 | ClassImp ( AliCdfJetFinder ) | |
50 | ||
51 | //______________________________________________________________________________ | |
52 | AliCdfJetFinder::AliCdfJetFinder(): | |
53 | AliJetFinder(), | |
54 | fHistos(0) | |
55 | { /* Constructor */ } | |
56 | ||
57 | //______________________________________________________________________________ | |
58 | AliCdfJetFinder::~AliCdfJetFinder() | |
59 | ||
60 | { | |
61 | // destructor | |
62 | cout << "Calling the destructor ... " << endl ; | |
63 | Reset(); | |
64 | cout << "Destructor called!" << endl; | |
65 | } | |
66 | ||
67 | //______________________________________________________________________________ | |
68 | void AliCdfJetFinder::CreateOutputObjects(TList *histos) | |
69 | { | |
70 | // Create the list of histograms. Only the list is owned. | |
71 | fHistos = histos; | |
72 | ||
73 | TH1F *h1 = new TH1F ("histo1", "no jets/event", 100, 0,20); | |
74 | h1->SetStats(kTRUE); | |
75 | h1->GetXaxis()->SetTitle("N_{jets}"); | |
76 | h1->GetYaxis()->SetTitle("#frac{dN}{dN_{jets}}"); | |
77 | h1->GetXaxis()->SetTitleColor(1); | |
78 | h1->SetMarkerStyle(kFullCircle); | |
79 | fHistos->Add(h1); | |
80 | ||
81 | TH1F *h2 = new TH1F ("histo2", "no part/leading jet", 40, 0,40); | |
82 | h2->SetStats(kTRUE); | |
83 | h2->GetXaxis()->SetTitle("N_{leading}"); | |
84 | h2->GetYaxis()->SetTitle("#frac{dN}{dN_{leading}}"); | |
85 | h2->GetXaxis()->SetTitleColor(1); | |
86 | h2->SetMarkerStyle(kFullCircle); | |
87 | fHistos->Add(h2); | |
88 | ||
89 | TProfile * h3 = new TProfile ("histo3","ProfileX of (pt,npart) of leading jets", 25, 0. ,25. , 0.,25. ) ; | |
90 | h3->SetStats(kTRUE); | |
91 | h3->GetXaxis()->SetTitle("P_{t}"); | |
92 | h3->GetYaxis()->SetTitle("N_{leading}"); | |
93 | h3->GetXaxis()->SetTitleColor(1); | |
94 | h3->SetMarkerStyle(kFullCircle); | |
95 | fHistos->Add(h3); | |
96 | ||
97 | TH1F *h4 = new TH1F ("histo4", "Charge momentum distribution for leading jet", 500, 0 , 5); | |
98 | fHistos->Add(h4); | |
99 | ||
100 | TProfile *h5 = new TProfile ("histo5", "ProfileX of N_{charge} vs dphi leading jet", 100 , 0. , 200. , 0 , 100 ); | |
101 | fHistos->Add(h5); | |
102 | ||
103 | TH1F *h6 = new TH1F ("histo6", " \"Transverse\" Pt Distribution ", 70, 0 , 14); | |
104 | fHistos->Add(h6); | |
105 | ||
106 | ||
107 | //cout << "CreateOutputObjetcs done!" << endl ; | |
108 | } | |
109 | ||
110 | //______________________________________________________________________________ | |
111 | void AliCdfJetFinder::FindJets() | |
112 | { | |
113 | //cout << "Begining of FindJets ..." << endl ; | |
114 | //1) Fill 1 array with momentum and initialise another array for indexes | |
115 | //2) Run algorithm | |
116 | // 2.1) Sort momentum array | |
117 | // 2.2) loop over arrays of TLorentzVectors | |
118 | // 2.3) flag as a possible jet | |
119 | //3) fill AliAODJet list | |
120 | Bool_t debug = kFALSE; | |
121 | ||
122 | TRefArray *refs = 0; | |
123 | Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader"); | |
124 | if (fromAod) { refs = fReader->GetReferences(); } | |
125 | ||
126 | AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader; | |
127 | Double_t radius = header->GetRadius(); // get Radius from header | |
128 | // cout << "Radius is : " << radius << endl ; | |
129 | ||
130 | TClonesArray * vectArray = fReader->GetMomentumArray() ; | |
131 | if ( vectArray == 0 ) { cout << "Could not get the momentum array" << endl; return; } | |
132 | ||
133 | Int_t nPart = vectArray->GetEntries() ; // n particles in this event; | |
134 | if ( nPart == 0 ) { if (debug) cout << "entries = 0 ; Event empty !!!" << endl ; return; } | |
135 | ||
136 | Int_t nJets = 0; // n jets in this event | |
137 | ||
138 | fJets->SetNinput ( nPart ) ; // number of input objects | |
139 | ||
140 | varContainer **vectParticle = new varContainer* [nPart]; // container for Particles | |
141 | varContainer **vectJet = new varContainer* [nPart]; // container for Jets | |
142 | ||
143 | Double_t *ptArray = new Double_t [nPart] ; // momentum array | |
144 | Int_t *idxArray = new Int_t [nPart] ; // index array of sorted pts | |
145 | ||
146 | // initialisation of momentum and index arrays | |
147 | // cout << "Filling idxArray && momArray for sorting " << endl; | |
148 | for ( Int_t i = 0 ; i < nPart ; i++ ) | |
149 | {// SORTING STEP :: ptArray with data from TClonesArray of TLorentzVector | |
150 | TLorentzVector * lv = (TLorentzVector*) vectArray->At(i); | |
151 | // INITIALISATION of local arrays for temporary storage | |
152 | varContainer *aParticle = new varContainer; | |
153 | aParticle->pt = lv->Pt(); | |
154 | aParticle->eta = lv->Eta(); | |
155 | aParticle->phi = lv->Phi(); | |
156 | aParticle->njet = -1; | |
157 | vectParticle[i] = aParticle; | |
158 | ||
159 | // initializing arrays | |
160 | idxArray [i] = 0 ; | |
161 | ptArray [i] = aParticle->pt ; | |
162 | } | |
163 | ||
164 | TMath::Sort ( nPart, ptArray, idxArray ) ; // get a sorted array of indexes with size = size of TClonesArray | |
165 | ||
166 | Double_t pt_jet = 0. , | |
167 | eta_jet = 0. , | |
168 | phi_jet = 0. , | |
169 | sum_eta = 0. , | |
170 | sum_phi = 0. ; // jet variables | |
171 | ||
172 | Int_t idxPtSort = 0 , // index of array of sorted pt indexes | |
173 | npartJet = 0 ; // number of particles in curent jet | |
174 | ||
175 | TBits lkup_table ( nPart ) ; // bit container of size npart | |
176 | ||
177 | while ( lkup_table.CountBits(0) != (UInt_t)nPart ) | |
178 | { // loop over particles in event until all flags are set | |
179 | UInt_t first_non_flagged = lkup_table.FirstNullBit(0) ; // set the index to the first NON flagged bit ; less conditions | |
180 | ||
181 | npartJet = 0 ; // reseting number of particles in jet | |
182 | pt_jet = 0.; | |
183 | eta_jet = 0.; | |
184 | phi_jet = 0.; | |
185 | sum_eta = 0.; | |
186 | sum_phi = 0.; | |
187 | ||
188 | for ( UInt_t i_part = first_non_flagged ; i_part < (UInt_t)nPart ; i_part++ ) | |
189 | {// iteration over particles in event | |
190 | // reseting variables | |
191 | ||
192 | // the loop is done over sorted array of pt | |
193 | idxPtSort = idxArray[i_part] ; // index entry of TLorentzVector from vect_arr | |
194 | ||
195 | if ( lkup_table.TestBitNumber(i_part) == 1 ) { continue; } // if 4vector is already flagged skip it | |
196 | ||
197 | //taking info from vectParticle ; | |
198 | Double_t pt_tmp = vectParticle[idxPtSort]->pt ; | |
199 | Double_t eta_tmp = vectParticle[idxPtSort]->eta ; | |
200 | Double_t phi_tmp = vectParticle[idxPtSort]->phi ; | |
201 | if (debug) printf(" particle %d: pt=%g\n", i_part, pt_tmp); | |
202 | ||
203 | // all angles are expressed in rad | |
204 | ||
205 | if ( TMath::Abs(eta_tmp) > 0.9 ) | |
206 | { | |
207 | lkup_table.SetBitNumber ( i_part, 1 ) ; // mark particle as used to be skipped | |
208 | continue ; | |
209 | } | |
210 | ||
211 | ||
212 | if ( i_part == first_non_flagged ) | |
213 | {// this is first particle in event; leading particle | |
214 | // initialise the jet variables with leading particle numbers | |
215 | ||
216 | npartJet++; // incrementing counter of particles in jet (nparJet == 1) (leading particle) | |
217 | pt_jet = pt_tmp ; | |
218 | eta_jet = eta_tmp ; | |
219 | phi_jet = phi_tmp ; | |
220 | if (debug) printf(" first part in jet: npartjet=%d pt_jet=%g\n", npartJet, pt_jet); | |
221 | ||
222 | lkup_table.SetBitNumber ( i_part, 1 ) ; // flag the 4vector | |
223 | ||
224 | vectParticle[idxPtSort]->njet = nJets ; // associate particle with current jet | |
225 | ||
226 | continue ; // skip to next particle | |
227 | } | |
228 | ||
229 | // condition to be in jet | |
230 | Double_t deta = eta_jet - eta_tmp ; | |
231 | Double_t dphi = phi_jet - phi_tmp ; | |
232 | ||
233 | if ( dphi < -TMath::Pi() ) { dphi = -dphi - 2.0 * TMath::Pi() ; } | |
234 | if ( dphi > TMath::Pi() ) { dphi = -dphi + 2.0 * TMath::Pi() ; } | |
235 | ||
236 | Double_t r_computed = TMath::Sqrt ( deta * deta + dphi * dphi ) ; | |
237 | ||
238 | Bool_t in_jet = ( r_computed <= radius ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0 | |
239 | ||
240 | if ( in_jet ) | |
241 | { // calculus of jet variables | |
242 | npartJet++; // incrementing counter of particles in jet | |
243 | ||
244 | pt_jet += pt_tmp ; | |
245 | sum_eta += (pt_tmp * eta_tmp) ; | |
246 | sum_phi += (pt_tmp * phi_tmp) ; | |
247 | eta_jet = sum_eta / pt_jet ; | |
248 | phi_jet = sum_phi / pt_jet ; | |
249 | ||
250 | lkup_table.SetBitNumber ( i_part, 1 ) ; // flag the 4vector | |
251 | vectParticle[idxPtSort]->njet = nJets ; // setting in particle list the associated jet | |
252 | if (debug) printf(" particle added to jet: npartjet=%d pt_jet=%g\n", npartJet, pt_jet); | |
253 | continue ; // skip to next particle | |
254 | } | |
255 | ||
256 | } | |
257 | // end of iteration over event; one jet definition | |
258 | ||
259 | varContainer *aJet = new varContainer; | |
260 | aJet->pt = pt_jet; | |
261 | aJet->eta = eta_jet; | |
262 | aJet->phi = phi_jet; | |
263 | aJet->njet = npartJet; | |
264 | vectJet[nJets++] = aJet; // store the jet (from jet 0 to nJets-1 ) | |
265 | if (debug) printf("=== current jet: npartjet=%d ptjet=%g\n", npartJet,pt_jet); | |
266 | ||
267 | // writing aod information | |
268 | Double_t px = 0., py = 0., pz = 0., en = 0.; // convert to 4-vector | |
269 | px = pt_jet * TMath::Cos ( phi_jet ) ; | |
270 | py = pt_jet * TMath::Sin ( phi_jet ) ; | |
271 | pz = pt_jet / TMath::Tan ( 2.0 * TMath::ATan ( TMath::Exp ( -eta_jet ) ) ) ; | |
272 | en = TMath::Sqrt ( px * px + py * py + pz * pz ); | |
273 | ||
274 | if (npartJet<2) continue; // do not add jets with less than 2 particles | |
275 | if (debug) cout << "Jet 4vect : " << "px = " << px << " ; py = " << py << " ; pz = " << pz << " ; E = " << en << endl; | |
276 | ||
277 | AliAODJet jet (px, py, pz, en); | |
278 | // cout << "Printing jet " << endl; | |
279 | if (debug) jet.Print(""); | |
280 | ||
281 | // cout << "Adding jet ... " ; | |
282 | AddJet(jet); | |
283 | // cout << "added \n" << endl; | |
284 | ||
285 | // if (fromAod) | |
286 | // { | |
287 | // for (Int_t parts = 0; parts < nPart; parts++ ) | |
288 | // { if (idx_jetT[parts] == nr_jet) {jet.AddTrack(refs->At(parts));} } | |
289 | // } | |
290 | ||
291 | ||
292 | } | |
293 | // end of while loop over particles ; ends when all particles were flagged as used and all jets defined | |
294 | ||
295 | ///////////////////////////////// | |
296 | //// END OF EVENT PARSING /// | |
297 | ///////////////////////////////// | |
298 | ||
299 | Int_t *jets_pt_idx = 0; // sorted array of jets pt | |
300 | Double_t *jets_pt = 0; // array of jets pts | |
301 | Int_t leading_jet_index = -1 ; // index of leading jet from vectJet | |
302 | Int_t part_leadjet = 0 ; // number of particles in leading jet | |
303 | Double_t pt_leadjet = 0. ; // pt of leading jet | |
304 | Double_t eta_leadjet = 0. ; // eta of leading jet | |
305 | Double_t phi_leadjet = 0. ; // phi of leading jet | |
306 | Int_t * idx_leadjet_part = 0; | |
307 | ||
308 | if (nJets > 0) | |
309 | { // if there is at least one jet in event | |
310 | jets_pt_idx = new Int_t [nJets] ; | |
311 | jets_pt = new Double_t [nJets] ; | |
312 | ||
313 | // filing the idx_ptjets array | |
314 | if (debug) printf("List of unsorted jets:\n"); | |
315 | for( Int_t i = 0 ; i < nJets ; i++ ) | |
316 | { | |
317 | jets_pt_idx [i] = 0 ; | |
318 | jets_pt [i] = vectJet[i]->pt ; | |
319 | if (debug) cout << " jet found: " << i << " npartjet=" << vectJet[i]->njet << " jets_pt[i]= " << jets_pt [i] << endl; | |
320 | } | |
321 | TMath::Sort ( nJets, jets_pt , jets_pt_idx ) ; // sorting pt of jets | |
322 | ||
323 | // selection of leading jet with nr of particles > 1 | |
324 | for( Int_t i = 0 ; i < nJets ; i++ ) | |
325 | { | |
326 | if ( vectJet[ jets_pt_idx[i] ]->njet > 1 ) | |
327 | { | |
328 | leading_jet_index = jets_pt_idx[i] ; | |
329 | part_leadjet = vectJet[ leading_jet_index ]->njet ; // number of particles in leading jet | |
330 | pt_leadjet = vectJet[ leading_jet_index ]->pt ; // pt of leading jet | |
331 | eta_leadjet = vectJet[ leading_jet_index ]->eta ; // eta of leading jet | |
332 | phi_leadjet = vectJet[ leading_jet_index ]->phi ; // phi of leading jet | |
333 | ||
334 | cout << "Leading jet: npart = " << part_leadjet | |
335 | << " ; pt = " << pt_leadjet | |
336 | << " ; phi = " << phi_leadjet | |
337 | << " ; eta = " << eta_leadjet | |
338 | << endl ; | |
339 | break ; | |
340 | } | |
341 | } | |
342 | // end of selection of leading jet | |
343 | ||
344 | // Filling an array with indexes of leading jet particles | |
345 | idx_leadjet_part = new Int_t [part_leadjet] ; | |
346 | Int_t counter = 0; | |
347 | if (debug) printf(" Searching particles with jet index %d\n", leading_jet_index); | |
348 | for( Int_t i = 0 ; i <nPart ; i++ ) | |
349 | { | |
350 | if ( vectParticle[i]->njet == leading_jet_index ) | |
351 | { | |
352 | idx_leadjet_part[counter++] = i ; | |
353 | if (debug) cout << " " << counter-1 << ": index=" << i << " pt=" << vectParticle[i]->pt << endl; | |
354 | } | |
355 | } | |
356 | if ( (counter-1) > part_leadjet ) { cout << " Counter > part_leadjet !!!!" << endl;} | |
357 | ||
358 | } | |
359 | // end of loop over DETECTED jets | |
360 | ||
361 | ||
362 | ||
363 | // Calculus of part distribution in leading jet | |
364 | Double_t z = 0. ; | |
365 | Double_t *z_part_ljet = new Double_t [ part_leadjet ] ; // array of z of particles in leading jet | |
366 | for( Int_t j = 0 ; j < part_leadjet ; j++ ) | |
367 | { | |
368 | Double_t z_j = vectParticle[idx_leadjet_part[j]]->pt ; | |
369 | z = z_j / pt_leadjet ; | |
370 | z_part_ljet [j] = z ; | |
371 | cout << "idx_leadjet_part[j] = " << idx_leadjet_part[j] | |
372 | << " p of particle = " << z_j | |
373 | << " pt lead jet = " << pt_leadjet | |
374 | << " Z = " << z << endl; | |
375 | } | |
376 | ||
377 | ||
378 | /* | |
379 | Double_t dphi_partLJ = 0. ; | |
380 | Double_t dphiArray [nPart]; | |
381 | for( Int_t part = 0 ; part < nPart ; part++ ) | |
382 | { | |
383 | dphi_partLJ = phi_leadjet - phi_partT [part] ; | |
384 | ||
385 | if ( dphi_partLJ < -TMath::Pi() ) { dphi_partLJ = -dphi_partLJ - 2.0 * TMath::Pi() ; } | |
386 | if ( dphi_partLJ > TMath::Pi() ) { dphi_partLJ = -dphi_partLJ + 2.0 * TMath::Pi() ; } | |
387 | ||
388 | dphi_partLJ = TMath::Abs(dphi_partLJ) ; | |
389 | ||
390 | dphiArray [part] = dphi_partLJ ; | |
391 | } | |
392 | */ | |
393 | // Fill nr_jet histogram | |
394 | if (fHistos) | |
395 | { | |
396 | // printf("FILLING histograms\n"); | |
397 | TH1F *h_jets = (TH1F*)fHistos->FindObject("histo1"); | |
398 | if (h_jets) h_jets->Fill(nJets); | |
399 | ||
400 | // There may be no leading jet with more that 2 particles ! | |
401 | TH1F *h_leadpart = (TH1F*)fHistos->FindObject("histo2"); | |
402 | if (part_leadjet && h_leadpart) h_leadpart->Fill(part_leadjet); | |
403 | ||
404 | TProfile * h_prof = (TProfile*)fHistos->FindObject("histo3"); | |
405 | if (part_leadjet && h_prof) h_prof->Fill(pt_leadjet,part_leadjet); | |
406 | ||
407 | TH1F *h_MD = (TH1F*)fHistos->FindObject("histo4"); | |
408 | for( Int_t k = 0 ; k < part_leadjet ; k++) | |
409 | { h_MD->Fill( z_part_ljet[k] ); } | |
410 | ||
411 | ||
412 | // TProfile *h_phi = (TProfile*)fHistos->At(4); | |
413 | // for( Int_t k = 0 ; k < nPart ; k++) | |
414 | // { h_phi->Fill( TMath::RadToDeg() * dphiArray [k] , nPart ) ; } | |
415 | ||
416 | ||
417 | // TH1F *h_tpd = (TH1F*)fHistos->At(5); | |
418 | // for( Int_t k = 0 ; k < nPart ; k++) | |
419 | // { | |
420 | // if ( (phi_partT[k] > (TMath::Pi()/3.)) && (phi_partT[k] < (2. * TMath::Pi()/3.)) ) | |
421 | // { h_tpd->Fill( pt_partT [k] ) ; } | |
422 | // } | |
423 | ||
424 | ||
425 | } | |
426 | ||
427 | ||
428 | // CLEANING SECTION | |
429 | for (Int_t i=0; i<nPart; i++) delete vectParticle[i]; | |
430 | delete [] vectParticle; | |
431 | for (Int_t i=0; i<nJets; i++) delete vectJet[i]; | |
432 | delete [] vectJet; | |
433 | delete [] ptArray; | |
434 | delete [] idxArray; | |
435 | ||
436 | if (z_part_ljet) delete [] z_part_ljet; | |
437 | if (idx_leadjet_part) delete [] idx_leadjet_part; | |
438 | if (jets_pt_idx) delete [] jets_pt_idx ; | |
439 | if (jets_pt) delete [] jets_pt ; | |
440 | ||
441 | vectArray->Delete() ; // TClonesArray of lorentz vectors | |
442 | } | |
443 | // end of FindJets | |
444 | ||
445 | ||
446 | //______________________________________________________________________________ | |
447 | void AliCdfJetFinder::FinishRun() | |
448 | {} | |
449 |