4f166ee35566d4c797f162bb92c45d2fe8a8f075
[u/mrichter/AliRoot.git] / JETAN / AliCdfJetFinder.cxx
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 usec, 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 // 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 /*
27 Changelog
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"
49 #include "AliAODJet.h"
50 #include "AliAODEvent.h"
51 #include "AliJetFinder.h"
52
53 #include "AliCdfJetFinder.h"
54 #include "AliCdfJetHeader.h"
55
56 ClassImp(AliCdfJetFinder)
57
58 //______________________________________________________________________________
59 AliCdfJetFinder::AliCdfJetFinder():
60     AliJetFinder(),
61     fHistos(0),
62     fDebug(0),
63     fFromAod(0),
64     fAODwrite(0),
65     fAODtracksWrite(0),
66     fAnalyseJets(0),
67     fRefArr (NULL),
68     fNJets(-9999),
69     fNPart(-9999),
70     fRadius(0.7),
71     fMinJetParticles(1),
72     fJetPtCut(0.),
73     fVectParticle(NULL),
74     fVectJet(NULL),
75     fPtArray(NULL),
76     fIdxArray(NULL)
77   {  /* Constructor */  }
78
79 //______________________________________________________________________________
80 AliCdfJetFinder::~AliCdfJetFinder()
81   {
82   // destructor
83   }
84
85 //______________________________________________________________________________
86 void 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 //______________________________________________________________________________
328 void 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".
338   if (fDebug) { printf("AliCDJetfinder::FindJets() %d \n", __LINE__ ); }
339 AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;
340
341   if (header)
342     {
343     fDebug            = header->IsDebugCDF();
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
353 if (fAODwrite)
354   {
355   fFromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
356   if (fFromAod) { fRefArr = fReader->GetReferences(); }
357   }
358
359 InitData();
360
361 if (!fNPart) { if (fDebug) {cout << "entries = 0 ; Event empty !!!" << endl ;} return; } // if event empty then exit
362
363 FindCones();
364
365 ComputeConesWeight();
366
367 if (fAODwrite) { 
368   if(fDebug)cout << "Writing AOD" << endl ; 
369   WriteJets();
370  }
371
372 if (fAnalyseJets) AnalizeJets();
373
374 Clean();
375
376 }
377
378 //______________________________________________________________________________
379 void AliCdfJetFinder::InitData()
380 {
381 // initialisation of variables and data members
382
383   TClonesArray * vectArray = fReader->GetMomentumArray() ;
384   if ( vectArray == 0 ) { cout << "Could not get the momentum array" << endl; return; }
385
386   fNPart = vectArray->GetEntries()  ; // n particles in this event;
387
388   if ( fNPart == 0 ) { return ; } // if event empty then exit
389
390   fVectParticle = new varContainer* [fNPart]; // container for Particles
391
392   fPtArray  = new Double_t [fNPart] ; // momentum array
393   fIdxArray = new Int_t    [fNPart] ; // index array of sorted pts
394
395   // initialisation of momentum and index arrays
396   for (  Int_t i = 0 ; i < fNPart ; i++ )
397     {// SORTING STEP :: fPtArray with data from TClonesArray of TLorentzVector
398     TLorentzVector * lv = (TLorentzVector*) vectArray->At(i);
399
400    // INITIALISATION of local arrays for temporary storage
401     varContainer *aParticle = new varContainer;
402     aParticle->pt   = lv->Pt();
403     aParticle->eta  = lv->Eta();
404     aParticle->phi  = TVector2::Phi_mpi_pi ( lv->Phi() ); // normalize to -pi,pi
405     aParticle->njet = -999;
406
407     fVectParticle[i] = aParticle;  // vector of Particles
408
409     // initializing arrays
410     fIdxArray [i] = -999 ;
411      fPtArray [i] = aParticle->pt ;
412     }
413
414   TMath::Sort ( fNPart, fPtArray, fIdxArray ) ; // get a sorted array of indexes with TClonesArray.Size()
415
416 }
417
418
419 //______________________________________________________________________________
420 void AliCdfJetFinder::FindCones()
421 {
422 // parsing of particles in event and estlabish jets (label them with jet index)
423
424   Double_t  ptSeed = 0. , etaSeed = 0. , phiSeed = 0. ; // leading particle params
425   Double_t pttmp = 0. , etatmp = 0. , phitmp = 0. ; // temporary variables to be used in various calculations
426   Double_t deta = 0. , dphi = 0. , dcomputed = 0. ;
427   Bool_t injet = 0 ;
428
429   fNJets = -1 ; // n jets in this event
430   Int_t idxPtSort = -1 ;  // index of array of sorted pt indexes
431
432   if (fDebug) { cout << "\n\n\n\n\n\n------------------\nBegin Event Analysis\n------------------\n\n" << endl ;}
433
434   if(fDebug)cout << "fNPart = " << fNPart << endl;
435
436   TBits lkupTable ( fNPart ) ;  // bit container ; 1-to-1 corespondence with fIdxArray
437
438   while ( lkupTable.CountBits() != (UInt_t)fNPart )
439     { // loop over particles in event until all flags are set
440     UInt_t firstnonflagged = lkupTable.FirstNullBit() ; // set the index to the first NON flagged bit ; less conditions
441
442     if(fDebug)cout << "\n\nfirst_non_flagged : " << firstnonflagged << endl;
443
444     ++fNJets; // incrementing the jet counter
445     if (fDebug) { printf("JET %d \n", fNJets); }
446
447     ptSeed = 0. ; etaSeed = 0. ; phiSeed = 0. ;  // reseting leading particle params
448
449     for (  UInt_t ipart = firstnonflagged ; ipart < (UInt_t)fNPart ; ipart++ )
450       {// iteration over particles in event
451       // the loop is done over sorted array of pt
452       idxPtSort = fIdxArray[ipart] ;  // index of particle ! fIdxArray is an index list pt sorted
453
454       if ( lkupTable.TestBitNumber(ipart) ) { continue; } // if 4vector is already flagged skip it
455
456       //init computed and used vars
457       pttmp = 0. ; etatmp = 0. ; phitmp = 0. ;
458       deta = 0. ; dphi = 0. ; dcomputed = 0. ; injet = 0 ;
459
460       //taking info from fVectParticle ;
461        pttmp = fVectParticle[idxPtSort]->pt ;
462       etatmp = fVectParticle[idxPtSort]->eta ;
463       phitmp = fVectParticle[idxPtSort]->phi ;
464
465       if ( ipart == firstnonflagged )
466         {// this is first particle in event; leading particle
467         // begin the search around this particle in a fRadius
468
469         // CENTRE OF THE JET
470         ptSeed = pttmp ; etaSeed = etatmp ; phiSeed = phitmp ; // seeding the jet with first particle idxPtSort
471
472         lkupTable.SetBitNumber ( ipart ) ; // flag the index of particle in lkup_table
473         fVectParticle[idxPtSort]->njet = fNJets ; // associate particle with current jet number
474
475         if (fDebug) { printf("\nLeading particle :: particle index = %d ;  at sorted index %d ; in jet %d \n", idxPtSort, ipart, fNJets); }
476         if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pttmp, etatmp, phitmp) ; }
477         if (fDebug) { lkupTable.Print() ;}
478
479         continue ; // skip to next particle
480         }
481
482       // condition to be in jet
483       deta = etatmp - etaSeed ;
484       dphi = TVector2::Phi_mpi_pi (phitmp - phiSeed) ; // computing dphi and normalizing to (0,2pi) interval in one step
485
486       dcomputed = TMath::Hypot(deta, dphi) ; // Distance(fRadius) to (eta,phi) seed
487
488       injet = ( ( fRadius - dcomputed ) >= 0.000000001 ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0
489
490       if ( injet )
491         { // calculus of jet variables
492         lkupTable.SetBitNumber ( ipart ) ;  // flag the index of particle in lkup_table
493         fVectParticle[idxPtSort]->njet = fNJets ; // setting in particle list the associated jet
494
495         if (fDebug) { printf("\njet particle :: particle index = %d ; at sorted index %d ; in jet %d ; found at radius %g ;  \n", idxPtSort, ipart, fNJets, dcomputed); }
496         if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pttmp, etatmp, phitmp) ; }
497         if (fDebug) { lkupTable.Print() ;}
498
499         continue ; // skip to next particle
500         }
501
502       }
503       // end of iteration over event; one jet definition of content ; jet parameters to be computed later
504     }
505 }
506
507
508 //______________________________________________________________________________
509 void AliCdfJetFinder::ComputeConesWeight()
510 {
511 // computing of jets Pt, Eta and Phi (centre of weight in (eta,phi) plane)
512 // rescan the vector of particles by identify them by asociate jet number for computing of weight centre
513
514 // JET CONTAINER
515 fVectJet      = new varContainer* [fNJets]; // container for Jets
516
517 Double_t ptJet, ptJet2 , etaJet , phiJet ; Int_t npartJet ;
518 Double_t pttmp = 0. , etatmp = 0. , phitmp = 0. ; // temporary variables to be used in various calculations
519 Int_t idxPtSort = -999 ;  // index of array of sorted pt indexes
520
521 for(  Int_t jet = 0 ; jet < fNJets ; jet++ )
522   {
523   if (fDebug) { printf("\n\n--- Computing weight of Jet %d \n", jet ); }
524   npartJet = 0 ; ptJet = 0. ; etaJet = 0. ; phiJet = 0. ; // reset variables for a new computation
525
526   for (  Int_t ipart = 0 ; ipart < fNPart ; ipart++ )
527     {// iteration over particles in event
528     // the loop is done over sorted array of pt
529     idxPtSort = fIdxArray[ipart] ;  // index of particle ! fIdxArray is an index list pt sorted
530
531     if ( fVectParticle[idxPtSort]->njet == jet )
532       {
533       ++npartJet; // incrementing the counter of jet particles
534
535       //taking info from fVectParticle ;
536        pttmp = fVectParticle[idxPtSort]->pt ;
537       etatmp = fVectParticle[idxPtSort]->eta ;
538       phitmp = TVector2::Phi_mpi_pi (fVectParticle[idxPtSort]->phi) ;
539
540 //      jet_new_angular_coordinate = jet_old_angular_coordinate * jet_old_pt / jet_new_pt +
541 //                                    part[i]_angular_coordinate * part[i]_pt/jet_new_pt
542
543       ptJet2 = ptJet + pttmp ;
544
545       etaJet = etaJet * ptJet / ptJet2 +  etatmp * pttmp / ptJet2 ;
546       phiJet = phiJet * ptJet / ptJet2 +  phitmp * pttmp / ptJet2 ;
547
548       ptJet = ptJet2 ;
549
550       }
551       // add a particle and recalculation of centroid
552     }
553     // end of 1 jet computation
554
555     varContainer *aJet = new varContainer;  // Jet container
556     aJet->pt = ptJet; aJet->eta = etaJet; aJet->phi = phiJet; aJet->njet = npartJet; // setting jet vars in container
557     fVectJet[jet] = aJet;   // store the number of the jet(fNJets) and increment afterwards
558
559     if (fDebug) { printf ("=== current jet %d : npartjet= %d ; pt_jet= %g ; eta_jet = %g ; phi_jet = %g \n\n\n",
560                                        jet,     npartJet,      ptJet,      etaJet,       phiJet ) ; }
561
562   }
563   //end loop over jets
564
565 }
566
567
568 //______________________________________________________________________________
569 void AliCdfJetFinder::WriteJets()  
570
571 // Writing AOD jets and AOD tracks
572
573 for(  Int_t jetnr = 0 ; jetnr < fNJets ; jetnr++ )
574   {
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[ jetnr ]->pt   ; // pt  of jet
578   eta = fVectJet[ jetnr ]->eta  ; // eta of jet
579   phi = fVectJet[ jetnr ]->phi  ; // phi of jet
580
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 );
585
586   AliAODJet jet (px, py, pz, en);
587
588
589   if (fDebug) jet.Print("");
590
591   if (fFromAod && fAODtracksWrite)
592     {
593       for (  Int_t jetTrack = 0; jetTrack < fNPart; jetTrack++ )
594         { if ( fVectParticle[jetTrack]->njet == jetnr ) { jet.AddTrack(fRefArr->At(jetTrack)) ; } }
595     }
596   // tracks REFs written in AOD
597   AddJet(jet);
598   }
599 //jets vector parsed and written to AOD
600 }
601
602
603 //______________________________________________________________________________
604 void AliCdfJetFinder::AnalizeJets()
605 {
606 // analyzing of jets and filling of histograms
607
608     const Double_t kPI = TMath::Pi();
609     
610   //persistent pointer to histo20
611   TH1F *hR = (TH1F*)fHistos->FindObject("histo20");
612
613   Int_t   *jetsptidx = 0;     // sorted array of jets pt
614   Double_t    *jetspt = 0;     // array of jets pts
615   Int_t leadingjetindex = -1 ;   // index of leading jet from fVectJet
616   Int_t partleadjet = 0 ; // number of particles in leading jet
617   Double_t   ptleadjet = 0. ; // pt  of leading jet
618   Double_t  etaleadjet = 0. ; // eta of leading jet
619   Double_t  phileadjet = 0. ; // phi of leading jet
620
621   jetsptidx = new Int_t    [fNJets] ;
622   jetspt    = new Double_t [fNJets] ;
623
624 //________________________________________________________________________________________
625 //  Jet sorting and finding the leading jet that coresponds to cuts in pt and multiplicity
626 //________________________________________________________________________________________
627
628   // filing the idx_ptjets array
629   if (fDebug) printf("List of unsorted jets:\n");
630   for(  Int_t i = 0 ; i < fNJets ; i++ )
631     {
632     jetsptidx [i] = 0 ;
633     jetspt    [i] = fVectJet[i]->pt ;
634     if (fDebug) { cout << "   jet found: " << i << " npartjet=" << fVectJet[i]->njet << " ; jets_pt = " << jetspt[i] << endl; }
635     }
636
637   TMath::Sort ( fNJets, jetspt , jetsptidx ) ; // sorting pt of jets
638
639   // selection of leading jet
640   // looping over jets searching for __first__ one that coresponds to cuts
641   for(  Int_t i = 0 ; i < fNJets ; i++ )
642     {
643     if ( ( fVectJet[ jetsptidx[i] ]->njet >= fMinJetParticles ) && ( fVectJet[ jetsptidx[i] ]->pt >= fJetPtCut ) )
644       {
645       leadingjetindex = jetsptidx[i] ;
646       partleadjet = fVectJet[ leadingjetindex ]->njet ; // number of particles in leading jet
647         ptleadjet = fVectJet[ leadingjetindex ]->pt   ; // pt  of leading jet
648        etaleadjet = fVectJet[ leadingjetindex ]->eta  ; // eta of leading jet
649        phileadjet = fVectJet[ leadingjetindex ]->phi  ; // phi of leading jet
650
651       if (fDebug)
652       { printf("Leading jet %d : npart= %d ; pt= %g ; eta = %g ; phi = %g \n", leadingjetindex, partleadjet, ptleadjet, etaleadjet, phileadjet ); }
653
654       break ;
655       }
656     }
657     // end of selection of leading jet
658
659
660
661 //////////////////////////////////////////////////
662 ////  Computing of values used in histograms
663 //////////////////////////////////////////////////
664
665 //___________________________________________________________________________
666 // pt_sum of all particles in event
667 //___________________________________________________________________________
668 cout << "Computing sum of pt in event" << endl ;
669 Double_t ptsumevent = 0.;
670 for (  Int_t i = 0 ; i< fNPart ; i++ ) { ptsumevent += fVectParticle[i]->pt ; }
671 printf ("Sum of all Pt in event : pt_sum_event = %g", ptsumevent) ;
672
673 //___________________________________________________________________________
674 // Filling an array with indexes of leading jet particles
675 //___________________________________________________________________________
676 Int_t * idxpartLJ = new Int_t [partleadjet] ;
677 Int_t counterpartleadjet = 0;
678
679 cout << "Filling an array with indexes of leading jet particles" << endl;
680
681 for( Int_t i = 0 ; i < fNPart ; i++ )
682   {
683   if ( fVectParticle[i]->njet == leadingjetindex )
684     {  idxpartLJ[counterpartleadjet++] = i ; }
685   }
686
687 if ( (counterpartleadjet-1) > partleadjet ) { cout << " Counter_part_lead_jet > part_leadjet !!!!" << endl;}
688
689
690 //___________________________________________________________________________
691 // Calculus of part distribution in leading jet
692 //___________________________________________________________________________
693 Double_t z = 0. ;
694 Double_t *zpartljet = new Double_t [ partleadjet ] ; // array of z of particles in leading jet
695
696 cout << "Entering loop of calculus of part distribution in leading jet" << endl ;
697
698 for( Int_t j = 0 ; j < partleadjet ; j++ )
699   {
700   Double_t zj = fVectParticle[idxpartLJ[j]]->pt ;
701   z =  zj / ptleadjet ;
702   zpartljet [j] = z ;
703   cout << "idx_leadjet_part[j] = " << idxpartLJ[j]
704       << " p of particle = " << zj
705       << " pt lead jet = " << ptleadjet
706       << " Z = " << z << endl;
707   }
708
709
710 //___________________________________________________________________________
711 // array of delta phi's between phi of particles and leading jet phi
712 //___________________________________________________________________________
713 cout << "array of delta phi's between phi of particles and leading jet phi" << endl;
714 Double_t dphipartLJ = 0. ;
715 Double_t *dphipartljet = new Double_t [fNPart];
716 for(  Int_t part = 0 ; part < fNPart ; part++ )
717   {
718   dphipartLJ = fVectParticle[part]->phi - phileadjet ;
719   dphipartLJ = TVector2::Phi_mpi_pi (dphipartLJ) ; // restrict the delta phi to (-pi,pi) interval
720   dphipartljet [part] = dphipartLJ ;
721   printf("part= %d ; dphi_partLJ = %g  \n", part, dphipartLJ );
722   }
723
724
725 //______________________________________________________________________________
726 //  Pt distribution for all particles
727 //______________________________________________________________________________
728 TH1F * hpt = (TH1F*)fHistos->FindObject("histo11");
729 if ( hpt ) { for (  Int_t i = 0 ; i < fNPart ; i++ ) { hpt->Fill( fVectParticle[i]->pt ); } }
730
731 //___________________________________________________________________________
732 // Recomputing of radius of particles in leading jet
733 //___________________________________________________________________________
734 if (fDebug) { printf("   Searching particles with jet index %d\n", leadingjetindex); }
735
736 Double_t ddeta = 0. , ddphi = 0. , rpart = 0. ;
737
738 for( Int_t j = 0 ; j < partleadjet ; j++ )
739   {
740   ddeta = etaleadjet - fVectParticle[idxpartLJ[j]]->eta;
741
742   Double_t phitmp = fVectParticle[idxpartLJ[j]]->phi ;
743
744   ddphi = TVector2::Phi_mpi_pi ( phileadjet - phitmp ) ; // restrict the delta phi to (-pi,pi) interval
745
746   rpart = TMath::Hypot (ddeta, ddphi) ;
747
748   printf ("Particle %d with Re-Computed radius = %f ", idxpartLJ[j], rpart) ;
749   if ( (rpart - fRadius) >= 0.00000001 )
750     { printf ("    bigger than selected radius of %f\n", fRadius ); }
751   else
752     { printf ("\n") ; }
753
754   if (hR) hR->Fill(rpart);
755
756   }
757
758
759
760 //_______________________________________________________________________
761 // Computing of radius that contain 80% of Leading Jet ( PT and multiplicity )
762 //_______________________________________________________________________
763 Double_t corepartleadjet = 0.8 * partleadjet ;
764 Double_t coreptleadjet = 0.8 * ptleadjet ;
765 Int_t countercorepart = 0 ;
766 Double_t countercorept = 0. ;
767 Int_t sortedindex = -1 ;
768
769 TProfile * hprof24 = (TProfile*)fHistos->FindObject("histo24");
770 TProfile * hprof25 = (TProfile*)fHistos->FindObject("histo25");
771
772 TProfile * hprof26 = (TProfile*)fHistos->FindObject("histo26");
773 TProfile * hprof27 = (TProfile*)fHistos->FindObject("histo27");
774 TProfile * hprof28 = (TProfile*)fHistos->FindObject("histo28");
775 TProfile * hprof29 = (TProfile*)fHistos->FindObject("histo29");
776
777
778 if ((hprof24) && (hprof25) && (hprof26) && (hprof27) && (hprof28) && (hprof29) )
779 {
780 for(  Int_t part = 0 ; part < fNPart ; part++ )
781   {
782   Double_t pttmp = 0. ; Double_t etatmp = 0. ; Double_t phitmp = 0. ; // temporary variables
783   Double_t dpart = 0. ;
784   sortedindex = fIdxArray[part] ;
785
786   if ( fVectParticle [ sortedindex ]->njet == leadingjetindex )
787     {
788      pttmp = fVectParticle[sortedindex]->pt ;
789     etatmp = fVectParticle[sortedindex]->eta ;
790     phitmp = fVectParticle[sortedindex]->phi ;
791
792     ++countercorepart ;
793     countercorept += pttmp ;
794
795     dpart = TMath::Hypot ( etaleadjet - etatmp, TVector2::Phi_mpi_pi (phileadjet - phitmp) ) ;
796
797     if ( countercorepart <=  corepartleadjet ) { hprof24->Fill(ptleadjet, dpart); }
798     if ( countercorept <= coreptleadjet ) { hprof25->Fill(ptleadjet, dpart); }
799
800     if (ptleadjet >  5.) { hprof26->Fill(dpart, countercorepart); hprof28->Fill(dpart, countercorept); }
801     if (ptleadjet > 30.) { hprof27->Fill(dpart, countercorepart); hprof29->Fill(dpart, countercorept); }
802
803     }
804   }
805 }
806
807   TH1F *hjetpt = (TH1F*)fHistos->FindObject("histo1");
808   TH1F *hjeteta = (TH1F*)fHistos->FindObject("histo2");
809   TH1F *hjetphi = (TH1F*)fHistos->FindObject("histo3");
810   TH1F *hjetnjet = (TH1F*)fHistos->FindObject("histo4");
811
812   for(  Int_t jet = 0 ; jet < fNJets ; jet++ )
813     {
814     if (hjetpt)   hjetpt   ->Fill ( fVectJet[jet]->pt   ) ;
815     if (hjeteta)  hjeteta  ->Fill ( fVectJet[jet]->eta  ) ;
816     if (hjetphi)  hjetphi  ->Fill ( fVectJet[jet]->phi  ) ;
817     if (hjetnjet) hjetnjet ->Fill ( fVectJet[jet]->njet ) ;
818     }
819
820   TH1F *hjets = (TH1F*)fHistos->FindObject("histo5");
821   if (hjets) hjets->Fill(fNJets);
822
823   TH1F *hleadpart = (TH1F*)fHistos->FindObject("histo6");
824   if (hleadpart) hleadpart->Fill(partleadjet);
825
826   TProfile * hprof = (TProfile*)fHistos->FindObject("histo7");
827   if (hprof) hprof->Fill(ptleadjet,partleadjet);
828
829   TH1F *hMD = (TH1F*)fHistos->FindObject("histo8");
830    for(  Int_t k = 0  ; k < partleadjet ; k++)
831      { hMD->Fill( zpartljet[k] ); }
832
833   TProfile * hphi = (TProfile*)fHistos->FindObject("histo9");
834     for(  Int_t k = 0  ; k < partleadjet ; k++)
835         { hphi->Fill( TMath::RadToDeg() * dphipartljet [k] , fNPart ) ; }
836
837   TProfile * htpd = (TProfile*)fHistos->FindObject("histo10");
838     for(  Int_t k = 0  ; k < fNPart ; k++)
839         { htpd->Fill( TMath::RadToDeg() * dphipartljet [k] , ptsumevent ) ; }
840
841
842   TProfile * hprof1 = (TProfile*)fHistos->FindObject("histo21");
843   if (hprof1) hprof1->Fill(ptleadjet, fNPart);
844
845   TProfile * hprof2 = (TProfile*)fHistos->FindObject("histo22");
846   if (hprof2) hprof2->Fill(ptleadjet, ptsumevent);
847
848   TProfile * hprof1toward = (TProfile*)fHistos->FindObject("histo21_toward");
849   TProfile * hprof1transverse = (TProfile*)fHistos->FindObject("histo21_transverse");
850   TProfile * hprof1away = (TProfile*)fHistos->FindObject("histo21_away");
851   TProfile * hprof2toward = (TProfile*)fHistos->FindObject("histo22_toward");
852   TProfile * hprof2transverse = (TProfile*)fHistos->FindObject("histo22_transverse");
853   TProfile * hprof2away = (TProfile*)fHistos->FindObject("histo22_away");
854   TH1F * hpttoward = (TH1F*)fHistos->FindObject("histo23_toward");
855   TH1F * hpttransverse = (TH1F*)fHistos->FindObject("histo23_transverse");
856   TH1F * hptaway = (TH1F*)fHistos->FindObject("histo23_away");
857
858   if ( (hprof1toward) && (hprof1transverse) && (hprof1away) && (hprof2toward) && (hprof2transverse) && (hprof2away) )
859     {
860     for( Int_t part = 0  ; part < fNPart ; part++)
861       {
862       Double_t ptpart = fVectParticle[part]->pt ; // pt of particle
863       if ( ( dphipartljet[part] >=0.) && ( dphipartljet[part] < kPI/3. ) )
864         {
865         hprof1toward->Fill( ptleadjet, fNPart );
866         hprof2toward->Fill( ptleadjet, ptsumevent);
867         hpttoward->Fill( ptpart );
868         }
869       else
870       if ( ( dphipartljet[part] >= (kPI/3.)) && ( dphipartljet[part] < (2.*kPI/3.)) )
871         {
872         hprof1transverse->Fill( ptleadjet, fNPart );
873         hprof2transverse->Fill( ptleadjet, ptsumevent);
874         hpttransverse->Fill( ptpart );
875         }
876       else
877       if ( ( dphipartljet[part] >= ( 2.*kPI/3.)) && ( dphipartljet[part] < kPI ) )
878         {
879         hprof1away->Fill( ptleadjet, fNPart );
880         hprof2away->Fill( ptleadjet, ptsumevent);
881         hptaway->Fill( ptpart );
882         }
883       }
884     }
885
886 }
887
888
889 //______________________________________________________________________________
890 void AliCdfJetFinder::Clean()
891 {
892 // CLEANING SECTION
893   delete [] fVectParticle;
894   delete [] fVectJet;
895   delete [] fPtArray;
896   delete [] fIdxArray;
897
898   Reset();
899 }
900
901 //______________________________________________________________________________
902 void AliCdfJetFinder::FinishRun()
903 {
904 // do i need this?
905 }
906