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 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 // 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 <vector>
  34 #include <iostream>
  35 #include <Riostream.h>
  36 #include <TFile.h>
  37 #include <TCanvas.h>
  38 #include <TROOT.h>
  39 #include <TClonesArray.h>
  40 #include <TLorentzVector.h>
  41 #include <TH1F.h>
  42 #include <TH2F.h>
  43 #include <TBits.h>
  44 #include <TArrayF.h>
  45 #include "AliCdfJetFinder.h"
  46 #include "AliCdfJetHeader.h"
  47 #include "AliJetReader.h"
  48 #include "AliJetReaderHeader.h"
  49 #include "AliJet.h"
  50 #include "AliAODJet.h"
  51 #include "AliAODEvent.h"
  52 #include "TProfile.h"
  53
  54
  55 ClassImp ( AliCdfJetFinder )
  56
  57 //______________________________________________________________________________
  58 AliCdfJetFinder::AliCdfJetFinder():
  59     AliJetFinder(),
  60     fHistos(0),
  61     fDebug(0),
  62     fFromAod(0),
  63     fAODwrite(0),
  64     fAODtracksWrite(0),
  65     fRefArr (NULL),
  66     fNJets(-9999),
  67     fNPart(-9999),
  68     fRadius(0.7),
  69     fMinJetParticles(1),
  70     fJetPtCut(0.),
  71     fVectParticle(NULL),
  72     fVectJet(NULL),
  73     fPtArray(NULL),
  74     fIdxArray(NULL)
  75   {  /* Constructor */  }
  76
  77 //______________________________________________________________________________
  78 AliCdfJetFinder::~AliCdfJetFinder()
  79
  80   {
  81   // destructor
  82   }
  83
  84 //______________________________________________________________________________
  85 void AliCdfJetFinder::CreateOutputObjects(TList *histos)
  86 {
  87   // Create the list of histograms. Only the list is owned.
  88   fHistos = histos;
  89
  90 //  gStyle->SetOptStat(11111111);
  91
  92   TH1F *h1 = new TH1F ("histo1", "Pt distribution of jets", 200, 0,200);  // 1GeV/bin
  93   h1->SetStats(kTRUE);
  94   h1->GetXaxis()->SetTitle("P_{T} of jets");
  95   h1->GetYaxis()->SetTitle("Number of jets");
  96   h1->GetXaxis()->SetTitleColor(1);
  97   h1->SetMarkerStyle(kFullCircle);
  98   fHistos->Add(h1);
  99
 100   TH1F *h2 = new TH1F ("histo2", "Eta distribution of jets", 240, -1.2,1.2); // 1 unit of rapidity / 100 bin
 101   h2->SetStats(kTRUE);
 102   h2->GetXaxis()->SetTitle("Eta of jets");
 103   h2->GetYaxis()->SetTitle("Number of jets");
 104   h2->GetXaxis()->SetTitleColor(1);
 105   h2->SetMarkerStyle(kFullCircle);
 106   fHistos->Add(h2);
 107
 108   TH1F *h3 = new TH1F ("histo3", "Phi distribution of jets", 400, -4,4);
 109   h3->SetStats(kTRUE);
 110   h3->GetXaxis()->SetTitle("Phi of jets");
 111   h3->GetYaxis()->SetTitle("Number of jets");
 112   h3->GetXaxis()->SetTitleColor(1);
 113   h3->SetMarkerStyle(kFullCircle);
 114   fHistos->Add(h3);
 115
 116   TH1F *h4 = new TH1F ("histo4", "Multiplicity of jets", 40, 0,40);  // 1 unit of multiplicity /bin
 117   h4->SetStats(kTRUE);
 118   h4->GetXaxis()->SetTitle("Particles in jets");
 119   h4->GetYaxis()->SetTitle("Number of jets");
 120   h4->GetXaxis()->SetTitleColor(1);
 121   h4->SetMarkerStyle(kFullCircle);
 122   fHistos->Add(h4);
 123
 124   TH1F *h5 = new TH1F ("histo5", "Distribution of jets in events", 100, 0,100);
 125   h5->SetStats(kTRUE);
 126   h5->GetXaxis()->SetTitle("Number of jets");
 127   h5->GetYaxis()->SetTitle("Number of events");
 128   h5->GetXaxis()->SetTitleColor(1);
 129   h5->SetMarkerStyle(kFullCircle);
 130   fHistos->Add(h5);
 131
 132   TH1F *h6 = new TH1F ("histo6", "Jet1 Charged Multiplicity Distribution", 30, 0,30);
 133   h6->SetStats(kTRUE);
 134   h6->GetXaxis()->SetTitle("N_{chg}");
 135   h6->GetYaxis()->SetTitle("Number of jets");
 136   h6->GetXaxis()->SetTitleColor(1);
 137   h6->SetMarkerStyle(kFullCircle);
 138   fHistos->Add(h6);
 139
 140   TProfile * h7 = new TProfile ("histo7","N_{chg}(jet1) vs P_{T}(charged jet1)", 200, 0. ,200. , 0.,200. ) ;
 141   h7->SetStats(kTRUE);
 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);
 146   fHistos->Add(h7);
 147
 148   TH1F *h8 = new TH1F ("histo8", "Charge momentum distribution for leading jet", 120, 0 , 1.2);
 149   h8->SetStats(kTRUE);
 150   h8->GetXaxis()->SetTitle("Jets");
 151   h8->GetYaxis()->SetTitle("Particle distribution");
 152   h8->GetXaxis()->SetTitleColor(1);
 153   h8->SetMarkerStyle(kFullCircle);
 154   fHistos->Add(h8);
 155
 156   TProfile *h9 = new TProfile ("histo9", "N_{chg} vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 20 );
 157   h9->SetStats(kTRUE);
 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);
 162   fHistos->Add(h9);
 163
 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);
 170   fHistos->Add(h10);
 171
 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);
 178   fHistos->Add(h11);
 179
 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);
 186   fHistos->Add(h20);
 187
 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);
 194   fHistos->Add(h21);
 195
 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);
 202   fHistos->Add(h22);
 203
 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);
 211
 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);
 219
 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);
 227
 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);
 235
 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);
 243
 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);
 251
 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);
 259
 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);
 267
 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);
 275
 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);
 282   fHistos->Add(h24);
 283
 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);
 290   fHistos->Add(h25);
 291
 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);
 298   fHistos->Add(h26);
 299
 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);
 306   fHistos->Add(h27);
 307
 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);
 314   fHistos->Add(h28);
 315
 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);
 322   fHistos->Add(h29);
 323
 324 }
 325
 326
 327 //______________________________________________________________________________
 328 /*
 329 Jet Algorithm:
 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".
 337
 338 */
 339 //______________________________________________________________________________
 340
 341 //______________________________________________________________________________
 342 void AliCdfJetFinder::FindJets()
 343 {
 344 AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;
 345
 346   if (header)
 347     {
 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
 354     }
 355   else
 356     { cout << "Header not found" << endl; return; }
 357
 358 // temporary until the other problems are resolved
 359 fAODwrite = 0 ;
 360
 361 if (fAODwrite)
 362   {
 363   fFromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
 364   if (fFromAod) { fRefArr = fReader->GetReferences(); }
 365   }
 366
 367 fFromAod = 0 ; // disable for the moment ; only ESD reading for now
 368
 369 InitData();
 370
 371 if (!fNPart) { if (fDebug) {cout << "entries = 0 ; Event empty !!!" << endl ;} return; } // if event empty then exit
 372
 373 FindCones();
 374
 375 ComputeConesWeight();
 376
 377 if (fAODwrite) { cout << "Writing AOD" << endl ; WriteJets(); }
 378
 379 AnalizeJets();
 380
 381 Clean();
 382
 383 }
 384
 385 //______________________________________________________________________________
 386 void AliCdfJetFinder::InitData()
 387 {
 388
 389   TClonesArray * vectArray = fReader->GetMomentumArray() ;
 390     if ( vectArray == 0 ) { cout << "Could not get the momentum array" << endl; return; }
 391
 392   fNPart = vectArray->GetEntries()  ; // n particles in this event;
 393
 394   if ( fNPart == 0 ) { return ; } // if event empty then exit
 395
 396   fJets->SetNinput ( fNPart ) ; // number of input objects
 397
 398   fVectParticle = new varContainer* [fNPart]; // container for Particles
 399   fVectJet      = new varContainer* [fNPart]; // container for Jets
 400
 401   fPtArray  = new Double_t [fNPart] ; // momentum array
 402   fIdxArray = new Int_t    [fNPart] ; // index array of sorted pts
 403
 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);
 408
 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;
 415
 416     fVectParticle[i] = aParticle;  // vector of Particles
 417
 418     // initializing arrays
 419     fIdxArray [i] = -999 ;
 420      fPtArray [i] = aParticle->pt ;
 421     }
 422
 423   TMath::Sort ( fNPart, fPtArray, fIdxArray ) ; // get a sorted array of indexes with TClonesArray.Size()
 424
 425 }
 426
 427
 428 //______________________________________________________________________________
 429 void AliCdfJetFinder::FindCones()
 430 {
 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. ;
 434   Bool_t in_jet = 0 ;
 435
 436   fNJets = -1 ; // n jets in this event
 437   Int_t idxPtSort = -1 ;  // index of array of sorted pt indexes
 438
 439   if (fDebug) { cout << "\n\n\n\n\n\n------------------\nBegin Event Analysis\n------------------\n\n" << endl ;}
 440
 441   cout << "fNPart = " << fNPart << endl;
 442
 443   TBits lkup_table ( fNPart ) ;  // bit container ; 1-to-1 corespondence with fIdxArray
 444
 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
 448
 449     cout << "\n\nfirst_non_flagged : " << first_non_flagged << endl;
 450
 451     ++fNJets; // incrementing the jet counter
 452     if (fDebug) { printf("JET %d \n", fNJets); }
 453
 454     pt_seed = 0. ; eta_seed = 0. ; phi_seed = 0. ;  // reseting leading particle params
 455
 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
 460
 461       if ( lkup_table.TestBitNumber(i_part) ) { continue; } // if 4vector is already flagged skip it
 462
 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 ;
 466
 467       //taking info from fVectParticle ;
 468        pt_tmp = fVectParticle[idxPtSort]->pt ;
 469       eta_tmp = fVectParticle[idxPtSort]->eta ;
 470       phi_tmp = fVectParticle[idxPtSort]->phi ;
 471
 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
 475
 476         // CENTRE OF THE JET
 477         pt_seed = pt_tmp ; eta_seed = eta_tmp ; phi_seed = phi_tmp ; // seeding the jet with first particle idxPtSort
 478
 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
 481
 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() ;}
 485
 486         continue ; // skip to next particle
 487         }
 488
 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
 492
 493       d_computed = Distance (deta, dphi) ; // Distance(fRadius) to (eta,phi) seed
 494
 495       in_jet = ( ( fRadius - d_computed ) >= 0.000000001 ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0
 496
 497       if ( in_jet )
 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
 501
 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() ;}
 505
 506         continue ; // skip to next particle
 507         }
 508
 509       }
 510       // end of iteration over event; one jet definition of content ; jet parameters to be computed later
 511     }
 512 }
 513
 514
 515 //______________________________________________________________________________
 516 void AliCdfJetFinder::ComputeConesWeight()
 517 {
 518 /*
 519       CALCULUS OF JETS Pt, Eta and Phi (centre of weight in (eta,phi) plane)
 520 */
 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
 523
 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
 527
 528 for(  Int_t jet = 0 ; jet < fNJets ; jet++ )
 529   {
 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
 532
 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
 537
 538     if ( fVectParticle[idxPtSort]->njet == jet )
 539       {
 540       ++npartJet; // incrementing the counter of jet particles
 541
 542       //taking info from fVectParticle ;
 543        pt_tmp = fVectParticle[idxPtSort]->pt ;
 544       eta_tmp = fVectParticle[idxPtSort]->eta ;
 545       phi_tmp = DeltaPhiNorm ( fVectParticle[idxPtSort]->phi ) ;
 546
 547       pt_jet += pt_tmp ;
 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) ;
 550
 551       }
 552       // add a particle and recalculation of centroid
 553     }
 554     // end of 1 jet computation
 555
 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
 559
 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 ) ; }
 562
 563   }
 564   //end loop over jets
 565
 566 }
 567
 568
 569 //______________________________________________________________________________
 570 void AliCdfJetFinder::WriteJets()
 571 { // Writing AOD jets and AOD tracks
 572
 573 /*  for(  Int_t jet_nr = 0 ; jet_nr < fNJets ; jet_nr++ )
 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[ jet_nr ]->pt   ; // pt  of jet
 578     eta = fVectJet[ jet_nr ]->eta  ; // eta of jet
 579     phi = fVectJet[ jet_nr ]->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     AddJet(jet);
 588
 589     if (fDebug) jet.Print("");
 590
 591     if (fromAod && AODtracksWrite)
 592       {
 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)) ; } }
 595       }
 596       // tracks REFs written in AOD
 597
 598     }*/
 599
 600 //jets vector parsed and written to AOD
 601
 602 }
 603
 604
 605 //______________________________________________________________________________
 606 void AliCdfJetFinder::AnalizeJets()
 607   {
 608
 609   //persistent pointer to histo20
 610   TH1F *h_r = (TH1F*)fHistos->FindObject("histo20");
 611
 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
 619
 620   jets_pt_idx = new Int_t    [fNJets] ;
 621   jets_pt     = new Double_t [fNJets] ;
 622
 623 //________________________________________________________________________________________
 624 //  Jet sorting and finding the leading jet that coresponds to cuts in pt and multiplicity
 625 //________________________________________________________________________________________
 626
 627   // filing the idx_ptjets array
 628   if (fDebug) printf("List of unsorted jets:\n");
 629   for(  Int_t i = 0 ; i < fNJets ; i++ )
 630     {
 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; }
 634     }
 635
 636   TMath::Sort ( fNJets, jets_pt , jets_pt_idx ) ; // sorting pt of jets
 637
 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++ )
 641     {
 642     if ( ( fVectJet[ jets_pt_idx[i] ]->njet >= fMinJetParticles ) && ( fVectJet[ jets_pt_idx[i] ]->pt >= fJetPtCut ) )
 643       {
 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
 649
 650       if (fDebug)
 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 ); }
 652
 653       break ;
 654       }
 655     }
 656     // end of selection of leading jet
 657
 658
 659
 660 //////////////////////////////////////////////////
 661 ////  Computing of values used in histograms
 662 //////////////////////////////////////////////////
 663
 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) ;
 671
 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;
 677
 678 cout << "Filling an array with indexes of leading jet particles" << endl;
 679
 680 for( Int_t i = 0 ; i < fNPart ; i++ )
 681   {
 682   if ( fVectParticle[i]->njet == leading_jet_index )
 683     {  idx_leadjet_part[counter_part_lead_jet++] = i ; }
 684   }
 685
 686 if ( (counter_part_lead_jet-1) > part_leadjet ) { cout << " Counter_part_lead_jet > part_leadjet !!!!" << endl;}
 687
 688
 689 //___________________________________________________________________________
 690 // Calculus of part distribution in leading jet
 691 //___________________________________________________________________________
 692 Double_t z = 0. ;
 693 Double_t *z_part_ljet = new Double_t [ part_leadjet ] ; // array of z of particles in leading jet
 694
 695 cout << "Entering loop of calculus of part distribution in leading jet" << endl ;
 696
 697 for( Int_t j = 0 ; j < part_leadjet ; j++ )
 698   {
 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;
 706   }
 707
 708
 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++ )
 716   {
 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 );
 721   }
 722
 723
 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 ); } }
 729
 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); }
 734
 735 Double_t ddeta = 0. , ddphi = 0. , r_part = 0. ;
 736
 737 for( Int_t j = 0 ; j < part_leadjet ; j++ )
 738   {
 739   ddeta = eta_leadjet - fVectParticle[idx_leadjet_part[j]]->eta;
 740
 741   Double_t phi_tmp = fVectParticle[idx_leadjet_part[j]]->phi ;
 742   phi_tmp = DeltaPhiNorm (phi_tmp);
 743
 744   ddphi = DeltaPhiNorm ( phi_leadjet - phi_tmp ) ; // restrict the delta phi to (-pi,pi) interval
 745
 746   r_part = Distance (ddeta, ddphi) ;
 747
 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 ); }
 751   else
 752     { printf ("\n") ; }
 753
 754   if (h_r) h_r->Fill(r_part);
 755
 756   }
 757
 758
 759
 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 ;
 768
 769 TProfile * h_prof_24 = (TProfile*)fHistos->FindObject("histo24");
 770 TProfile * h_prof_25 = (TProfile*)fHistos->FindObject("histo25");
 771
 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");
 776
 777
 778 if ((h_prof_24) && (h_prof_25) && (h_prof_26) && (h_prof_27) && (h_prof_28) && (h_prof_29) )
 779 {
 780 for(  Int_t part = 0 ; part < fNPart ; part++ )
 781   {
 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] ;
 785
 786   if ( fVectParticle [ sorted_index ]->njet == leading_jet_index )
 787     {
 788      pt_tmp = fVectParticle[sorted_index]->pt ;
 789     eta_tmp = fVectParticle[sorted_index]->eta ;
 790     phi_tmp = fVectParticle[sorted_index]->phi ;
 791
 792     ++counter_core_part ;
 793     counter_core_pt += pt_tmp ;
 794
 795     d_part = Distance ( eta_leadjet - eta_tmp, phi_leadjet - phi_tmp ) ;
 796
 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); }
 799
 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); }
 802
 803     }
 804   }
 805 }
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 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");
 820
 821   for(  Int_t jet = 0 ; jet < fNJets ; jet++ )
 822     {
 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 ) ;
 827     }
 828
 829   TH1F *h_jets = (TH1F*)fHistos->FindObject("histo5");
 830   if (h_jets) h_jets->Fill(fNJets);
 831
 832   TH1F *h_leadpart = (TH1F*)fHistos->FindObject("histo6");
 833   if (h_leadpart) h_leadpart->Fill(part_leadjet);
 834
 835   TProfile * h_prof = (TProfile*)fHistos->FindObject("histo7");
 836   if (h_prof) h_prof->Fill(pt_leadjet,part_leadjet);
 837
 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] ); }
 841
 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 ) ; }
 845
 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 ) ; }
 849
 850
 851   TProfile * h_prof1 = (TProfile*)fHistos->FindObject("histo21");
 852   if (h_prof1) h_prof1->Fill(pt_leadjet, fNPart);
 853
 854   TProfile * h_prof2 = (TProfile*)fHistos->FindObject("histo22");
 855   if (h_prof2) h_prof2->Fill(pt_leadjet, pt_sum_event);
 856
 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");
 866
 867
 868
 869   if ( (h_prof1_toward) && (h_prof1_transverse) && (h_prof1_away) && (h_prof2_toward) && (h_prof2_transverse) && (h_prof2_away) )
 870     {
 871     for( Int_t part = 0  ; part < fNPart ; part++)
 872       {
 873       Double_t pt_part = fVectParticle[part]->pt ; // pt of particle
 874       if ( ( dphi_part_ljet[part] >=0.) && ( dphi_part_ljet[part] < pi/3. ) )
 875         {
 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 );
 879         }
 880       else
 881       if ( ( dphi_part_ljet[part] >= (pi/3.)) && ( dphi_part_ljet[part] < (2.*pi/3.)) )
 882         {
 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 );
 886         }
 887       else
 888       if ( ( dphi_part_ljet[part] >= ( 2.*pi/3.)) && ( dphi_part_ljet[part] < pi ) )
 889         {
 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 );
 893         }
 894       }
 895     }
 896
 897
 898
 899
 900
 901
 902
 903 }
 904
 905
 906 //______________________________________________________________________________
 907 void AliCdfJetFinder::Clean()
 908 {
 909
 910
 911   // CLEANING SECTION
 912   delete [] fVectParticle;
 913   delete [] fVectJet;
 914   delete [] fPtArray;
 915   delete [] fIdxArray;
 916
 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 ;
 921
 922 //  vectArray->Delete() ; // TClonesArray of lorentz vectors
 923
 924   Reset();
 925
 926
 927
 928 }
 929
 930
 931
 932
 933 //______________________________________________________________________________
 934 void AliCdfJetFinder::FinishRun()
 935 {}
 936
 937