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