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 /* $Id$ */
  17
  18 //---------------------------------------------------------------------
  19 // Jet Finder based on CDF algortihm
  20 // Charged jet evolution and the underlying event in proton-antiproton collisions at 1.8 TeV
  21 // Physical Review D, vol. 65, Issue 9, id. 092002
  22 // http://www.phys.ufl.edu/~rfield/cdf/chgjet/chgjet_intro.html
  23 // Authors : Adrian.Sevcenco@cern.ch (adriansev@spacescience.ro )
  24 //           Daniel.Felea@cern.ch    (dfelea@spacescience.ro)
  25 //           Ciprian.Mihai.Mitu@cern.ch  (mcm@spacescience.ro)
  26 //           magali.estienne@subatech.in2p3.fr &
  27 //           alexandre.shabetai@cern.ch (Modification of the input object (reader/finder splitting))
  28 // ** 2011
  29 // Modified accordingly to reader/finder splitting and new handling of neutral information
  30 //---------------------------------------------------------------------
  31
  32 #include <Riostream.h>
  33 #include <TMath.h>
  34 #include <TBits.h>
  35 #include <TFile.h>
  36 #include <TH1F.h>
  37 #include <TProfile.h>
  38 #include <TVector2.h>
  39
  40 #include "AliAODJet.h"
  41 #include "AliJetFinder.h"
  42 #include "AliJetCalTrk.h"
  43 #include "AliCdfJetFinder.h"
  44 #include "AliCdfJetHeader.h"
  45
  46 ClassImp(AliCdfJetFinder)
  47
  48 //______________________________________________________________________________
  49 AliCdfJetFinder::AliCdfJetFinder():
  50   AliJetFinder(),
  51   fHistos(0),
  52   fAODwrite(0),
  53   fAODtracksWrite(0),
  54   fAnalyseJets(0),
  55   fNJets(0),
  56   fNPart(0),
  57   fNInC(0),
  58   fNInN(0),
  59   fRadius(0.7),
  60   fMinJetParticles(1),
  61   fJetPtCut(0.),
  62   fVectParticle(NULL),
  63   fVectJet(NULL),
  64   fPtArray(NULL),
  65   fIdxArray(NULL)
  66 {
  67   // Default constructor
  68 }
  69
  70 //______________________________________________________________________________
  71 AliCdfJetFinder::~AliCdfJetFinder()
  72 {
  73   // destructor
  74   Clean();
  75
  76 }
  77
  78 //______________________________________________________________________________
  79 void AliCdfJetFinder::CreateOutputObjects(TList * const histos)
  80 {
  81   // Create the list of histograms. Only the list is owned.
  82   fHistos = histos;
  83
  84   //  gStyle->SetOptStat(11111111);
  85
  86   TH1F *h1 = new TH1F ("histo1", "Pt distribution of jets", 200, 0,200);  // 1GeV/bin
  87   h1->SetStats(kTRUE);
  88   h1->GetXaxis()->SetTitle("P_{T} of jets");
  89   h1->GetYaxis()->SetTitle("Number of jets");
  90   h1->GetXaxis()->SetTitleColor(1);
  91   h1->SetMarkerStyle(kFullCircle);
  92   fHistos->Add(h1);
  93
  94   TH1F *h2 = new TH1F ("histo2", "Eta distribution of jets", 240, -1.2,1.2); // 1 unit of rapidity / 100 bin
  95   h2->SetStats(kTRUE);
  96   h2->GetXaxis()->SetTitle("Eta of jets");
  97   h2->GetYaxis()->SetTitle("Number of jets");
  98   h2->GetXaxis()->SetTitleColor(1);
  99   h2->SetMarkerStyle(kFullCircle);
 100   fHistos->Add(h2);
 101
 102   TH1F *h3 = new TH1F ("histo3", "Phi distribution of jets", 400, -4,4);
 103   h3->SetStats(kTRUE);
 104   h3->GetXaxis()->SetTitle("Phi of jets");
 105   h3->GetYaxis()->SetTitle("Number of jets");
 106   h3->GetXaxis()->SetTitleColor(1);
 107   h3->SetMarkerStyle(kFullCircle);
 108   fHistos->Add(h3);
 109
 110   TH1F *h4 = new TH1F ("histo4", "Multiplicity of jets", 40, 0,40);  // 1 unit of multiplicity /bin
 111   h4->SetStats(kTRUE);
 112   h4->GetXaxis()->SetTitle("Particles in jets");
 113   h4->GetYaxis()->SetTitle("Number of jets");
 114   h4->GetXaxis()->SetTitleColor(1);
 115   h4->SetMarkerStyle(kFullCircle);
 116   fHistos->Add(h4);
 117
 118   TH1F *h5 = new TH1F ("histo5", "Distribution of jets in events", 100, 0,100);
 119   h5->SetStats(kTRUE);
 120   h5->GetXaxis()->SetTitle("Number of jets");
 121   h5->GetYaxis()->SetTitle("Number of events");
 122   h5->GetXaxis()->SetTitleColor(1);
 123   h5->SetMarkerStyle(kFullCircle);
 124   fHistos->Add(h5);
 125
 126   TH1F *h6 = new TH1F ("histo6", "Jet1 Charged Multiplicity Distribution", 30, 0,30);
 127   h6->SetStats(kTRUE);
 128   h6->GetXaxis()->SetTitle("N_{chg}");
 129   h6->GetYaxis()->SetTitle("Number of jets");
 130   h6->GetXaxis()->SetTitleColor(1);
 131   h6->SetMarkerStyle(kFullCircle);
 132   fHistos->Add(h6);
 133
 134   TProfile * h7 = new TProfile ("histo7","N_{chg}(jet1) vs P_{T}(charged jet1)", 200, 0. ,200. , 0.,200. ) ;
 135   h7->SetStats(kTRUE);
 136   h7->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 137   h7->GetYaxis()->SetTitle("<N_{chg}(jet1)> in 1 GeV/c bin");
 138   h7->GetXaxis()->SetTitleColor(1);
 139   h7->SetMarkerStyle(kFullCircle);
 140   fHistos->Add(h7);
 141
 142   TH1F *h8 = new TH1F ("histo8", "Charge momentum distribution for leading jet", 120, 0 , 1.2);
 143   h8->SetStats(kTRUE);
 144   h8->GetXaxis()->SetTitle("Jets");
 145   h8->GetYaxis()->SetTitle("Particle distribution");
 146   h8->GetXaxis()->SetTitleColor(1);
 147   h8->SetMarkerStyle(kFullCircle);
 148   fHistos->Add(h8);
 149
 150   TProfile *h9 = new TProfile ("histo9", "N_{chg} vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 20 );
 151   h9->SetStats(kTRUE);
 152   h9->GetXaxis()->SetTitle("#Delta#phi (degrees)");
 153   h9->GetYaxis()->SetTitle("<N_{chg}> in 3.6 degree bin");
 154   h9->GetXaxis()->SetTitleColor(1);
 155   h9->SetMarkerStyle(kFullCircle);
 156   fHistos->Add(h9);
 157
 158   TProfile *h10 = new TProfile ("histo10", "P_{T} sum vs the Azimuthal Angle from Charged Jet1", 50 , 0. , 180. , 0 , 100 );
 159   h10->SetStats(kTRUE);
 160   h10->GetXaxis()->SetTitle("#Delta#phi (degrees)");
 161   h10->GetYaxis()->SetTitle("<P_{T} sum> in 3.6 degree bin");
 162   h10->GetXaxis()->SetTitleColor(1);
 163   h10->SetMarkerStyle(kFullCircle);
 164   fHistos->Add(h10);
 165
 166   TH1F *h11 = new TH1F ("histo11", " \"Transverse\" Pt Distribution ", 70, 0 , 14);
 167   h11->SetStats(kTRUE);
 168   h11->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 169   h11->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
 170   h11->GetXaxis()->SetTitleColor(1);
 171   h11->SetMarkerStyle(kFullCircle);
 172   fHistos->Add(h11);
 173
 174   TH1F *h20 = new TH1F ("histo20", "Distribution of R in leading jet", 400, 0.,4.);
 175   h20->SetStats(kTRUE);
 176   h20->GetXaxis()->SetTitle("R [formula]");
 177   h20->GetYaxis()->SetTitle("dN/dR");
 178   h20->GetXaxis()->SetTitleColor(1);
 179   h20->SetMarkerStyle(kFullCircle);
 180   fHistos->Add(h20);
 181
 182   TProfile * h21 = new TProfile ("histo21","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 30. ) ;
 183   h21->SetStats(kTRUE);
 184   h21->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 185   h21->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
 186   h21->GetXaxis()->SetTitleColor(1);
 187   h21->SetMarkerStyle(kFullCircle);
 188   fHistos->Add(h21);
 189
 190   TProfile * h22 = new TProfile ("histo22","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
 191   h22->SetStats(kTRUE);
 192   h22->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 193   h22->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
 194   h22->GetXaxis()->SetTitleColor(1);
 195   h22->SetMarkerStyle(kFullCircle);
 196   fHistos->Add(h22);
 197
 198   TProfile * h21Toward = new TProfile ("histo21_toward","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
 199   h21Toward->SetStats(kTRUE);
 200   h21Toward->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 201   h21Toward->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
 202   h21Toward->GetXaxis()->SetTitleColor(1);
 203   h21Toward->SetMarkerStyle(kFullCircle);
 204   fHistos->Add(h21Toward);
 205
 206   TProfile * h21Transverse = new TProfile ("histo21_transverse","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
 207   h21Transverse->SetStats(kTRUE);
 208   h21Transverse->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 209   h21Transverse->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
 210   h21Transverse->GetXaxis()->SetTitleColor(1);
 211   h21Transverse->SetMarkerStyle(kFullCircle);
 212   fHistos->Add(h21Transverse);
 213
 214   TProfile * h21Away = new TProfile ("histo21_away","N_{chg}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0., 50. , 0., 12. ) ;
 215   h21Away->SetStats(kTRUE);
 216   h21Away->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 217   h21Away->GetYaxis()->SetTitle("<N_{chg}(in the event - including jet1)> in 1 GeV/c bin");
 218   h21Away->GetXaxis()->SetTitleColor(1);
 219   h21Away->SetMarkerStyle(kFullCircle);
 220   fHistos->Add(h21Away);
 221
 222   TProfile * h22Toward = new TProfile ("histo22_toward","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
 223   h22Toward->SetStats(kTRUE);
 224   h22Toward->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 225   h22Toward->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
 226   h22Toward->GetXaxis()->SetTitleColor(1);
 227   h22Toward->SetMarkerStyle(kFullCircle);
 228   fHistos->Add(h22Toward);
 229
 230   TProfile * h22Transverse = new TProfile ("histo22_transverse","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
 231   h22Transverse->SetStats(kTRUE);
 232   h22Transverse->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 233   h22Transverse->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
 234   h22Transverse->GetXaxis()->SetTitleColor(1);
 235   h22Transverse->SetMarkerStyle(kFullCircle);
 236   fHistos->Add(h22Transverse);
 237
 238   TProfile * h22Away = new TProfile ("histo22_away","PT_{sum}(in the event - including jet1) vs P_{T}(charged jet1)", 200, 0. , 50. , 0., 50. ) ;
 239   h22Away->SetStats(kTRUE);
 240   h22Away->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 241   h22Away->GetYaxis()->SetTitle("<PT_{sum}(in the event - including jet1)> in 1 GeV/c bin");
 242   h22Away->GetXaxis()->SetTitleColor(1);
 243   h22Away->SetMarkerStyle(kFullCircle);
 244   fHistos->Add(h22Away);
 245
 246   TH1F *h23Toward = new TH1F ("histo23_toward","'Toward' Pt Distribution of charged particles", 200, 0., 14.);
 247   h23Toward->SetStats(kTRUE);
 248   h23Toward->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
 249   h23Toward->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
 250   h23Toward->GetXaxis()->SetTitleColor(1);
 251   h23Toward->SetMarkerStyle(kFullCircle);
 252   fHistos->Add(h23Toward);
 253
 254   TH1F *h23Transverse = new TH1F ("histo23_transverse","'Transverse' Pt Distribution of charged particles", 200, 0., 14.);
 255   h23Transverse->SetStats(kTRUE);
 256   h23Transverse->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
 257   h23Transverse->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
 258   h23Transverse->GetXaxis()->SetTitleColor(1);
 259   h23Transverse->SetMarkerStyle(kFullCircle);
 260   fHistos->Add(h23Transverse);
 261
 262   TH1F *h23Away = new TH1F ("histo23_away","'Away' Pt Distribution of charged particles", 200, 0., 14.);
 263   h23Away->SetStats(kTRUE);
 264   h23Away->GetXaxis()->SetTitle("P_{T} (charged) (GeV/c)");
 265   h23Away->GetYaxis()->SetTitle("dN_{chg}/dP_{T} (1/GeV/c)");
 266   h23Away->GetXaxis()->SetTitleColor(1);
 267   h23Away->SetMarkerStyle(kFullCircle);
 268   fHistos->Add(h23Away);
 269
 270   TProfile * h24 = new TProfile ("histo24","Jet1 Size vs P_{T}(charged jet1)", 200, 0., 50. , 0., 0.5) ;
 271   h24->SetStats(kTRUE);
 272   h24->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 273   h24->GetYaxis()->SetTitle("<R(chgjet1)> in 1 GeV/c bin");
 274   h24->GetXaxis()->SetTitleColor(1);
 275   h24->SetMarkerStyle(kFullCircle);
 276   fHistos->Add(h24);
 277
 278   TProfile * h25 = new TProfile ("histo25","Jet1 Size vs P_{T}(charged jet1)", 200, 0., 50. , 0., 0.5) ;
 279   h25->SetStats(kTRUE);
 280   h25->GetXaxis()->SetTitle("P_{T} (charged jet1) (GeV/c)");
 281   h25->GetYaxis()->SetTitle("<R(chgjet1)> in 1 GeV/c bin");
 282   h25->GetXaxis()->SetTitleColor(1);
 283   h25->SetMarkerStyle(kFullCircle);
 284   fHistos->Add(h25);
 285
 286   TProfile *h26 = new TProfile ("histo26", "N_{chg} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0., 0.8);
 287   h26->SetStats(kTRUE);
 288   h26->GetXaxis()->SetTitle("Distance R");
 289   h26->GetYaxis()->SetTitle("<N_{chg}> in 0.02 bin");
 290   h26->GetXaxis()->SetTitleColor(1);
 291   h26->SetMarkerStyle(kFullCircle);
 292   fHistos->Add(h26);
 293
 294   TProfile *h27 = new TProfile ("histo27", "N_{chg} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0., 0.8);
 295   h27->SetStats(kTRUE);
 296   h27->GetXaxis()->SetTitle("Distance R");
 297   h27->GetYaxis()->SetTitle("<N_{chg}> in 0.02 bin");
 298   h27->GetXaxis()->SetTitleColor(1);
 299   h27->SetMarkerStyle(kFullCircle);
 300   fHistos->Add(h27);
 301
 302   TProfile *h28 = new TProfile ("histo28", "PT_{sum} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0.01, 10.);
 303   h28->SetStats(kTRUE);
 304   h28->GetXaxis()->SetTitle("Distance R");
 305   h28->GetYaxis()->SetTitle("<PT_{sum} (GeV/c)> in 0.02 bin");
 306   h28->GetXaxis()->SetTitleColor(1);
 307   h28->SetMarkerStyle(kFullCircle);
 308   fHistos->Add(h28);
 309
 310   TProfile *h29 = new TProfile ("histo29", "PT_{sum} vs the Distance R from Charged Jet1", 30, 0., 0.6, 0.01, 10.);
 311   h29->SetStats(kTRUE);
 312   h29->GetXaxis()->SetTitle("Distance R");
 313   h29->GetYaxis()->SetTitle("<PT_{sum} (GeV/c)> in 0.02 bin");
 314   h29->GetXaxis()->SetTitleColor(1);
 315   h29->SetMarkerStyle(kFullCircle);
 316   fHistos->Add(h29);
 317
 318 }
 319
 320 //______________________________________________________________________________
 321 void AliCdfJetFinder::FindJets()
 322 {
 323   // Jet Algorithm:
 324   //  * Order all charged particles according to their PT.
 325   //  * Start with the highest PT particle and include in the "jet" all particles within the "radius" R = 0.7
 326   //     (considering each particle in the order of decreasing PT and recalculating the centroid of the jet after
 327   //     each new particle is added to the jet).
 328   //  * Go to the next highest PT particle (not already included in a jet) and include in the "jet" all particles
 329   //     (not already included in a jet) within the radius R =0.7.
 330   //  * Continue until all particles are in a "jet".
 331   if (fDebug) { printf("AliCDJetfinder::FindJets() %d \n", __LINE__ ); }
 332   AliCdfJetHeader *header = (AliCdfJetHeader*)fHeader;
 333
 334   if (header)
 335     {
 336       fDebug            = header->GetDebug();
 337       fAODwrite         = header->IsAODwrite() ;       // write jets to AOD
 338       fAODtracksWrite   = header->IsAODtracksWrite() ; // write jet tracks to AOD
 339       fRadius           = header->GetRadius();      // get Radius from jet finder header
 340       fMinJetParticles  = header->GetMinPartJet (); // get minimum multiplicity of an jet
 341       fJetPtCut         = header->GetJetPtCut ();   // get minimum of jet pt
 342       fAnalyseJets      = header->GetAnalyseJets(); // get analyse jet
 343     }
 344   else
 345     { cout << "Header not found" << endl; return; }
 346
 347   InitData();
 348
 349   if (!fNPart) {
 350     if (fDebug) {
 351       cout << "entries = 0 ; Event empty !!!" << endl ;
 352     }
 353     // no need to call clean, InitData does not
 354     // create pointers if npart == 0
 355     return;
 356   } // if event empty then exit
 357
 358   FindCones();
 359
 360   ComputeConesWeight();
 361
 362   if (fAODwrite) {
 363     if(fDebug)cout << "Writing AOD" << endl ;
 364     WriteJets();
 365   }
 366
 367   if (fAnalyseJets) AnalizeJets();
 368
 369   Clean();
 370
 371 }
 372
 373 //______________________________________________________________________________
 374 void AliCdfJetFinder::InitData()
 375 {
 376   // initialisation of variables and data members
 377
 378   if( fHeader->GetDebug() && fCalTrkEvent->GetNCalTrkTracks()  == 0) { cout << "No charged tracks found" << endl; }
 379
 380   fNPart = fCalTrkEvent->GetNCalTrkTracks() ;
 381
 382   if ( fCalTrkEvent->GetNCalTrkTracks() ) { return; } // if event empty then exit
 383
 384   fVectParticle = new varContainer* [fNPart]; // container for Particles
 385
 386   fPtArray  = new Double_t [fCalTrkEvent->GetNCalTrkTracks()] ;
 387   fIdxArray = new Int_t    [fCalTrkEvent->GetNCalTrkTracks()] ; // index array of sorted pts
 388
 389   // initialisation of momentum and index arrays
 390   for (  Int_t i = 0 ; i <fCalTrkEvent->GetNCalTrkTracks() ; i++ )
 391     {// SORTING STEP :: fPtArray with data from CalTrkTracks
 392
 393       // INITIALISATION of local arrays for temporary storage
 394       varContainer *aParticle = new varContainer;
 395       aParticle->pt   = fCalTrkEvent->GetCalTrkTrack(i)->GetPt();
 396       aParticle->eta  = fCalTrkEvent->GetCalTrkTrack(i)->GetEta();
 397       aParticle->phi  = TVector2::Phi_mpi_pi ( fCalTrkEvent->GetCalTrkTrack(i)->GetPhi() ); // normalize to -pi,pi
 398       aParticle->njet = -999;
 399
 400       fVectParticle[i] = aParticle;  // vector of Particles
 401
 402       // initializing arrays
 403       fIdxArray [i] = -999 ;
 404       fPtArray [i] = aParticle->pt ;
 405     }
 406
 407   TMath::Sort ( fNPart, fPtArray, fIdxArray ) ; // get a sorted array of indexes
 408
 409 }
 410
 411 //______________________________________________________________________________
 412 void AliCdfJetFinder::FindCones()
 413 {
 414   // parsing of particles in event and estlabish jets (label them with jet index)
 415
 416   Double_t  ptSeed = 0. , etaSeed = 0. , phiSeed = 0. ; // leading particle params
 417   Double_t pttmp = 0. , etatmp = 0. , phitmp = 0. ; // temporary variables to be used in various calculations
 418   Double_t deta = 0. , dphi = 0. , dcomputed = 0. ;
 419   Bool_t injet = 0 ;
 420
 421   fNJets = -1 ; // n jets in this event
 422   Int_t idxPtSort = -1 ;  // index of array of sorted pt indexes
 423
 424   if (fDebug) { cout << "\n\n\n\n\n\n------------------\nBegin Event Analysis\n------------------\n\n" << endl ;}
 425
 426   if(fDebug)cout << "fNPart = " << fNPart << endl;
 427
 428   TBits lkupTable ( fNPart ) ;  // bit container ; 1-to-1 corespondence with fIdxArray
 429
 430   while ( lkupTable.CountBits() != (UInt_t)fNPart )
 431     { // loop over particles in event until all flags are set
 432       UInt_t firstnonflagged = lkupTable.FirstNullBit() ; // set the index to the first NON flagged bit ; less conditions
 433
 434       if(fDebug)cout << "\n\nfirst_non_flagged : " << firstnonflagged << endl;
 435
 436       ++fNJets; // incrementing the jet counter
 437       if (fDebug) { printf("JET %d \n", fNJets); }
 438
 439       ptSeed = 0. ; etaSeed = 0. ; phiSeed = 0. ;  // reseting leading particle params
 440
 441       for (  UInt_t ipart = firstnonflagged ; ipart < (UInt_t)fNPart ; ipart++ )
 442         {// iteration over particles in event
 443           // the loop is done over sorted array of pt
 444           idxPtSort = fIdxArray[ipart] ;  // index of particle ! fIdxArray is an index list pt sorted
 445
 446           if ( lkupTable.TestBitNumber(ipart) ) { continue; } // if 4vector is already flagged skip it
 447
 448           //init computed and used vars
 449           pttmp = 0. ; etatmp = 0. ; phitmp = 0. ;
 450           deta = 0. ; dphi = 0. ; dcomputed = 0. ; injet = 0 ;
 451
 452           //taking info from fVectParticle ;
 453           pttmp = fVectParticle[idxPtSort]->pt ;
 454           etatmp = fVectParticle[idxPtSort]->eta ;
 455           phitmp = fVectParticle[idxPtSort]->phi ;
 456
 457           if ( ipart == firstnonflagged )
 458             {// this is first particle in event; leading particle
 459               // begin the search around this particle in a fRadius
 460
 461               // CENTRE OF THE JET
 462               ptSeed = pttmp ; etaSeed = etatmp ; phiSeed = phitmp ; // seeding the jet with first particle idxPtSort
 463
 464               lkupTable.SetBitNumber ( ipart ) ; // flag the index of particle in lkup_table
 465               fVectParticle[idxPtSort]->njet = fNJets ; // associate particle with current jet number
 466
 467               if (fDebug) { printf("\nLeading particle :: particle index = %d ;  at sorted index %d ; in jet %d \n", idxPtSort, ipart, fNJets); }
 468               if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pttmp, etatmp, phitmp) ; }
 469               if (fDebug) { lkupTable.Print() ;}
 470
 471               continue ; // skip to next particle
 472             }
 473
 474           // condition to be in jet
 475           deta = etatmp - etaSeed ;
 476           dphi = TVector2::Phi_mpi_pi (phitmp - phiSeed) ; // computing dphi and normalizing to (0,2pi) interval in one step
 477
 478           dcomputed = TMath::Hypot(deta, dphi) ; // Distance(fRadius) to (eta,phi) seed
 479
 480           injet = ( ( fRadius - dcomputed ) >= 0.000000001 ) ? 1 : 0 ; // if r_computed is within jet_r in_jet == 1 else 0
 481
 482           if ( injet )
 483             { // calculus of jet variables
 484               lkupTable.SetBitNumber ( ipart ) ;  // flag the index of particle in lkup_table
 485               fVectParticle[idxPtSort]->njet = fNJets ; // setting in particle list the associated jet
 486
 487               if (fDebug) { printf("\njet particle :: particle index = %d ; at sorted index %d ; in jet %d ; found at radius %g ;  \n", idxPtSort, ipart, fNJets, dcomputed); }
 488               if (fDebug) { printf("pt= %g ; eta= %g ; phi = %g \n", pttmp, etatmp, phitmp) ; }
 489               if (fDebug) { lkupTable.Print() ;}
 490
 491               continue ; // skip to next particle
 492             }
 493
 494         }
 495       // end of iteration over event; one jet definition of content ; jet parameters to be computed later
 496     }
 497
 498 }
 499
 500 //______________________________________________________________________________
 501 void AliCdfJetFinder::ComputeConesWeight()
 502 {
 503   // computing of jets Pt, Eta and Phi (centre of weight in (eta,phi) plane)
 504   // rescan the vector of particles by identify them by asociate jet number for computing of weight centre
 505
 506   // JET CONTAINER
 507   fVectJet      = new varContainer* [fNJets]; // container for Jets
 508
 509   Double_t ptJet, ptJet2 , etaJet , phiJet ; Int_t npartJet ;
 510   Double_t pttmp = 0. , etatmp = 0. , phitmp = 0. ; // temporary variables to be used in various calculations
 511   Int_t idxPtSort = -999 ;  // index of array of sorted pt indexes
 512
 513   for(  Int_t jet = 0 ; jet < fNJets ; jet++ )
 514     {
 515       if (fDebug) { printf("\n\n--- Computing weight of Jet %d \n", jet ); }
 516       npartJet = 0 ; ptJet = 0. ; etaJet = 0. ; phiJet = 0. ; // reset variables for a new computation
 517
 518       for (  Int_t ipart = 0 ; ipart < fNPart ; ipart++ )
 519         {// iteration over particles in event
 520           // the loop is done over sorted array of pt
 521           idxPtSort = fIdxArray[ipart] ;  // index of particle ! fIdxArray is an index list pt sorted
 522
 523           if ( fVectParticle[idxPtSort]->njet == jet )
 524             {
 525               ++npartJet; // incrementing the counter of jet particles
 526
 527               //taking info from fVectParticle ;
 528               pttmp = fVectParticle[idxPtSort]->pt ;
 529               etatmp = fVectParticle[idxPtSort]->eta ;
 530               phitmp = TVector2::Phi_mpi_pi (fVectParticle[idxPtSort]->phi) ;
 531
 532               //      jet_new_angular_coordinate = jet_old_angular_coordinate * jet_old_pt / jet_new_pt +
 533               //                                    part[i]_angular_coordinate * part[i]_pt/jet_new_pt
 534
 535               ptJet2 = ptJet + pttmp ;
 536
 537               etaJet = etaJet * ptJet / ptJet2 +  etatmp * pttmp / ptJet2 ;
 538               phiJet = phiJet * ptJet / ptJet2 +  phitmp * pttmp / ptJet2 ;
 539
 540               ptJet = ptJet2 ;
 541
 542             }
 543           // add a particle and recalculation of centroid
 544         }
 545       // end of 1 jet computation
 546
 547       varContainer *aJet = new varContainer;  // Jet container
 548       aJet->pt = ptJet; aJet->eta = etaJet; aJet->phi = phiJet; aJet->njet = npartJet; // setting jet vars in container
 549       fVectJet[jet] = aJet;   // store the number of the jet(fNJets) and increment afterwards
 550
 551       if (fDebug) { printf ("=== current jet %d : npartjet= %d ; pt_jet= %g ; eta_jet = %g ; phi_jet = %g \n\n\n",
 552                             jet,     npartJet,      ptJet,      etaJet,       phiJet ) ; }
 553
 554     }
 555   //end loop over jets
 556
 557 }
 558
 559 //______________________________________________________________________________
 560 void AliCdfJetFinder::WriteJets()
 561 {
 562   // Writing AOD jets and AOD tracks
 563
 564   for(  Int_t jetnr = 0 ; jetnr < fNJets ; jetnr++ )
 565     {
 566       Double_t pt = 0., eta = 0., phi = 0., // jet variables
 567         px = 0., py = 0., pz = 0., en = 0.; // convert to 4-vector
 568       pt  = fVectJet[ jetnr ]->pt   ; // pt  of jet
 569       eta = fVectJet[ jetnr ]->eta  ; // eta of jet
 570       phi = fVectJet[ jetnr ]->phi  ; // phi of jet
 571
 572       px = pt * TMath::Cos ( phi ) ;
 573       py = pt * TMath::Sin ( phi ) ;
 574       pz = pt / TMath::Tan ( 2.0 * TMath::ATan ( TMath::Exp ( -eta ) ) ) ;
 575       en = TMath::Sqrt ( px * px + py * py + pz * pz );
 576
 577       AliAODJet jet (px, py, pz, en);
 578
 579
 580       if (fDebug) jet.Print("");
 581
 582       if (fAODtracksWrite)
 583         {
 584           for (  Int_t jetTrack = 0; jetTrack < fCalTrkEvent->GetNCalTrkTracks(); jetTrack++ )
 585             {
 586               // The first if condition below has to be checked
 587               if ( fVectParticle[jetTrack]->njet == jetnr ) { jet.AddTrack(fCalTrkEvent->GetCalTrkTrack(jetTrack)->GetTrackObject()) ; }
 588             }
 589         }
 590       // tracks REFs written in AOD
 591       AddJet(jet);
 592
 593     }
 594   //jets vector parsed and written to AOD
 595 }
 596
 597 //______________________________________________________________________________
 598 void AliCdfJetFinder::AnalizeJets()
 599 {
 600   // analyzing of jets and filling of histograms
 601
 602   const Double_t kPI = TMath::Pi();
 603
 604   //persistent pointer to histo20
 605   TH1F *hR = (TH1F*)fHistos->FindObject("histo20");
 606
 607   Int_t   *jetsptidx = 0;     // sorted array of jets pt
 608   Double_t    *jetspt = 0;     // array of jets pts
 609   Int_t leadingjetindex = -1 ;   // index of leading jet from fVectJet
 610   Int_t partleadjet = 0 ; // number of particles in leading jet
 611   Double_t   ptleadjet = 0. ; // pt  of leading jet
 612   Double_t  etaleadjet = 0. ; // eta of leading jet
 613   Double_t  phileadjet = 0. ; // phi of leading jet
 614
 615   jetsptidx = new Int_t    [fNJets] ;
 616   jetspt    = new Double_t [fNJets] ;
 617
 618   //________________________________________________________________________________________
 619   //  Jet sorting and finding the leading jet that coresponds to cuts in pt and multiplicity
 620   //________________________________________________________________________________________
 621
 622   // filing the idx_ptjets array
 623   if (fDebug) printf("List of unsorted jets:\n");
 624   for(  Int_t i = 0 ; i < fNJets ; i++ )
 625     {
 626       jetsptidx [i] = 0 ;
 627       jetspt    [i] = fVectJet[i]->pt ;
 628       if (fDebug) { cout << "   jet found: " << i << " npartjet=" << fVectJet[i]->njet << " ; jets_pt = " << jetspt[i] << endl; }
 629     }
 630
 631   TMath::Sort ( fNJets, jetspt , jetsptidx ) ; // sorting pt of jets
 632
 633   // selection of leading jet
 634   // looping over jets searching for __first__ one that coresponds to cuts
 635   for(  Int_t i = 0 ; i < fNJets ; i++ )
 636     {
 637       if ( ( fVectJet[ jetsptidx[i] ]->njet >= fMinJetParticles ) && ( fVectJet[ jetsptidx[i] ]->pt >= fJetPtCut ) )
 638         {
 639           leadingjetindex = jetsptidx[i] ;
 640           partleadjet = fVectJet[ leadingjetindex ]->njet ; // number of particles in leading jet
 641           ptleadjet = fVectJet[ leadingjetindex ]->pt   ; // pt  of leading jet
 642           etaleadjet = fVectJet[ leadingjetindex ]->eta  ; // eta of leading jet
 643           phileadjet = fVectJet[ leadingjetindex ]->phi  ; // phi of leading jet
 644
 645           if (fDebug)
 646             { printf("Leading jet %d : npart= %d ; pt= %g ; eta = %g ; phi = %g \n", leadingjetindex, partleadjet, ptleadjet, etaleadjet, phileadjet ); }
 647
 648           break ;
 649         }
 650     }
 651   // end of selection of leading jet
 652
 653
 654
 655   //////////////////////////////////////////////////
 656   ////  Computing of values used in histograms
 657   //////////////////////////////////////////////////
 658
 659   //___________________________________________________________________________
 660   // pt_sum of all particles in event
 661   //___________________________________________________________________________
 662   if (fDebug) cout << "Computing sum of pt in event" << endl ;
 663   Double_t ptsumevent = 0.;
 664   for (  Int_t i = 0 ; i< fNPart ; i++ ) { ptsumevent += fVectParticle[i]->pt ; }
 665   if (fDebug) printf ("Sum of all Pt in event : pt_sum_event = %g", ptsumevent) ;
 666
 667   //___________________________________________________________________________
 668   // Filling an array with indexes of leading jet particles
 669   //___________________________________________________________________________
 670   Int_t * idxpartLJ = new Int_t [partleadjet] ;
 671   Int_t counterpartleadjet = 0;
 672
 673   if (fDebug) cout << "Filling an array with indexes of leading jet particles" << endl;
 674
 675   for( Int_t i = 0 ; i < fNPart ; i++ )
 676     {
 677       if ( fVectParticle[i]->njet == leadingjetindex )
 678         {  idxpartLJ[counterpartleadjet++] = i ; }
 679     }
 680
 681   if ( (counterpartleadjet-1) > partleadjet ) { cout << " Counter_part_lead_jet > part_leadjet !!!!" << endl;}
 682
 683
 684   //___________________________________________________________________________
 685   // Calculus of part distribution in leading jet
 686   //___________________________________________________________________________
 687   Double_t z = 0. ;
 688   Double_t *zpartljet = new Double_t [ partleadjet ] ; // array of z of particles in leading jet
 689
 690   if (fDebug) cout << "Entering loop of calculus of part distribution in leading jet" << endl ;
 691
 692   for( Int_t j = 0 ; j < partleadjet ; j++ )
 693     {
 694       Double_t zj = fVectParticle[idxpartLJ[j]]->pt ;
 695       z =  zj / ptleadjet ;
 696       zpartljet [j] = z ;
 697       if (fDebug) cout << "idx_leadjet_part[j] = " << idxpartLJ[j]
 698                        << " p of particle = " << zj
 699                        << " pt lead jet = " << ptleadjet
 700                        << " Z = " << z << endl;
 701     }
 702
 703
 704   //___________________________________________________________________________
 705   // array of delta phi's between phi of particles and leading jet phi
 706   //___________________________________________________________________________
 707   if (fDebug) cout << "array of delta phi's between phi of particles and leading jet phi" << endl;
 708   Double_t dphipartLJ = 0. ;
 709   Double_t *dphipartljet = new Double_t [fNPart];
 710   for(  Int_t part = 0 ; part < fNPart ; part++ )
 711     {
 712       dphipartLJ = fVectParticle[part]->phi - phileadjet ;
 713       dphipartLJ = TVector2::Phi_mpi_pi (dphipartLJ) ; // restrict the delta phi to (-pi,pi) interval
 714       dphipartljet [part] = dphipartLJ ;
 715       if (fDebug) printf("part= %d ; dphi_partLJ = %g  \n", part, dphipartLJ );
 716     }
 717
 718
 719   //______________________________________________________________________________
 720   //  Pt distribution for all particles
 721   //______________________________________________________________________________
 722   TH1F * hpt = (TH1F*)fHistos->FindObject("histo11");
 723   if ( hpt ) { for (  Int_t i = 0 ; i < fNPart ; i++ ) { hpt->Fill( fVectParticle[i]->pt ); } }
 724
 725   //___________________________________________________________________________
 726   // Recomputing of radius of particles in leading jet
 727   //___________________________________________________________________________
 728   if (fDebug) { printf("   Searching particles with jet index %d\n", leadingjetindex); }
 729
 730   Double_t ddeta = 0. , ddphi = 0. , rpart = 0. ;
 731
 732   for( Int_t j = 0 ; j < partleadjet ; j++ )
 733     {
 734       ddeta = etaleadjet - fVectParticle[idxpartLJ[j]]->eta;
 735
 736       Double_t phitmp = fVectParticle[idxpartLJ[j]]->phi ;
 737
 738       ddphi = TVector2::Phi_mpi_pi ( phileadjet - phitmp ) ; // restrict the delta phi to (-pi,pi) interval
 739
 740       rpart = TMath::Hypot (ddeta, ddphi) ;
 741
 742       if (fDebug) printf ("Particle %d with Re-Computed radius = %f ", idxpartLJ[j], rpart) ;
 743       if ( (rpart - fRadius) >= 0.00000001 )
 744         { if (fDebug) printf ("    bigger than selected radius of %f\n", fRadius ); }
 745       else
 746         { if (fDebug) printf ("\n") ; }
 747
 748       if (hR) hR->Fill(rpart);
 749
 750     }
 751
 752
 753
 754   //_______________________________________________________________________
 755   // Computing of radius that contain 80% of Leading Jet ( PT and multiplicity )
 756   //_______________________________________________________________________
 757   Double_t corepartleadjet = 0.8 * partleadjet ;
 758   Double_t coreptleadjet = 0.8 * ptleadjet ;
 759   Int_t countercorepart = 0 ;
 760   Double_t countercorept = 0. ;
 761   Int_t sortedindex = -1 ;
 762
 763   TProfile * hprof24 = (TProfile*)fHistos->FindObject("histo24");
 764   TProfile * hprof25 = (TProfile*)fHistos->FindObject("histo25");
 765   TProfile * hprof26 = (TProfile*)fHistos->FindObject("histo26");
 766   TProfile * hprof27 = (TProfile*)fHistos->FindObject("histo27");
 767   TProfile * hprof28 = (TProfile*)fHistos->FindObject("histo28");
 768   TProfile * hprof29 = (TProfile*)fHistos->FindObject("histo29");
 769
 770
 771   if ((hprof24) && (hprof25) && (hprof26) && (hprof27) && (hprof28) && (hprof29) )
 772     {
 773       for(  Int_t part = 0 ; part < fNPart ; part++ )
 774         {
 775           Double_t pttmp = 0. ; Double_t etatmp = 0. ; Double_t phitmp = 0. ; // temporary variables
 776           Double_t dpart = 0. ;
 777           sortedindex = fIdxArray[part] ;
 778
 779           if ( fVectParticle [ sortedindex ]->njet == leadingjetindex )
 780             {
 781               pttmp = fVectParticle[sortedindex]->pt ;
 782               etatmp = fVectParticle[sortedindex]->eta ;
 783               phitmp = fVectParticle[sortedindex]->phi ;
 784
 785               ++countercorepart ;
 786               countercorept += pttmp ;
 787
 788               dpart = TMath::Hypot ( etaleadjet - etatmp, TVector2::Phi_mpi_pi (phileadjet - phitmp) ) ;
 789
 790               if ( countercorepart <=  corepartleadjet ) { hprof24->Fill(ptleadjet, dpart); }
 791               if ( countercorept <= coreptleadjet ) { hprof25->Fill(ptleadjet, dpart); }
 792
 793               if (ptleadjet >  5.) { hprof26->Fill(dpart, countercorepart); hprof28->Fill(dpart, countercorept); }
 794               if (ptleadjet > 30.) { hprof27->Fill(dpart, countercorepart); hprof29->Fill(dpart, countercorept); }
 795
 796             }
 797         }
 798     }
 799
 800   TH1F *hjetpt = (TH1F*)fHistos->FindObject("histo1");
 801   TH1F *hjeteta = (TH1F*)fHistos->FindObject("histo2");
 802   TH1F *hjetphi = (TH1F*)fHistos->FindObject("histo3");
 803   TH1F *hjetnjet = (TH1F*)fHistos->FindObject("histo4");
 804
 805   for(  Int_t jet = 0 ; jet < fNJets ; jet++ )
 806     {
 807       if (hjetpt)   hjetpt   ->Fill ( fVectJet[jet]->pt   ) ;
 808       if (hjeteta)  hjeteta  ->Fill ( fVectJet[jet]->eta  ) ;
 809       if (hjetphi)  hjetphi  ->Fill ( fVectJet[jet]->phi  ) ;
 810       if (hjetnjet) hjetnjet ->Fill ( fVectJet[jet]->njet ) ;
 811     }
 812
 813   TH1F *hjets = (TH1F*)fHistos->FindObject("histo5");
 814   if (hjets) hjets->Fill(fNJets);
 815
 816   TH1F *hleadpart = (TH1F*)fHistos->FindObject("histo6");
 817   if (hleadpart) hleadpart->Fill(partleadjet);
 818
 819   TProfile * hprof = (TProfile*)fHistos->FindObject("histo7");
 820   if (hprof) hprof->Fill(ptleadjet,partleadjet);
 821
 822   TH1F *hMD = (TH1F*)fHistos->FindObject("histo8");
 823   for(  Int_t k = 0  ; k < partleadjet ; k++)
 824     { hMD->Fill( zpartljet[k] ); }
 825
 826   TProfile * hphi = (TProfile*)fHistos->FindObject("histo9");
 827   for(  Int_t k = 0  ; k < partleadjet ; k++)
 828     { hphi->Fill( TMath::RadToDeg() * dphipartljet [k] , fNPart ) ; }
 829
 830   TProfile * htpd = (TProfile*)fHistos->FindObject("histo10");
 831   for(  Int_t k = 0  ; k < fNPart ; k++)
 832     { htpd->Fill( TMath::RadToDeg() * dphipartljet [k] , ptsumevent ) ; }
 833
 834
 835   TProfile * hprof1 = (TProfile*)fHistos->FindObject("histo21");
 836   if (hprof1) hprof1->Fill(ptleadjet, fNPart);
 837
 838   TProfile * hprof2 = (TProfile*)fHistos->FindObject("histo22");
 839   if (hprof2) hprof2->Fill(ptleadjet, ptsumevent);
 840
 841   TProfile * hprof1toward = (TProfile*)fHistos->FindObject("histo21_toward");
 842   TProfile * hprof1transverse = (TProfile*)fHistos->FindObject("histo21_transverse");
 843   TProfile * hprof1away = (TProfile*)fHistos->FindObject("histo21_away");
 844   TProfile * hprof2toward = (TProfile*)fHistos->FindObject("histo22_toward");
 845   TProfile * hprof2transverse = (TProfile*)fHistos->FindObject("histo22_transverse");
 846   TProfile * hprof2away = (TProfile*)fHistos->FindObject("histo22_away");
 847   TH1F * hpttoward = (TH1F*)fHistos->FindObject("histo23_toward");
 848   TH1F * hpttransverse = (TH1F*)fHistos->FindObject("histo23_transverse");
 849   TH1F * hptaway = (TH1F*)fHistos->FindObject("histo23_away");
 850
 851   if ( (hprof1toward) && (hprof1transverse) && (hprof1away) && (hprof2toward) && (hprof2transverse) && (hprof2away) )
 852     {
 853       for( Int_t part = 0  ; part < fNPart ; part++)
 854         {
 855           Double_t ptpart = fVectParticle[part]->pt ; // pt of particle
 856           if ( ( dphipartljet[part] >=0.) && ( dphipartljet[part] < kPI/3. ) )
 857             {
 858               hprof1toward->Fill( ptleadjet, fNPart );
 859               hprof2toward->Fill( ptleadjet, ptsumevent);
 860               hpttoward->Fill( ptpart );
 861             }
 862           else
 863             if ( ( dphipartljet[part] >= (kPI/3.)) && ( dphipartljet[part] < (2.*kPI/3.)) )
 864               {
 865                 hprof1transverse->Fill( ptleadjet, fNPart );
 866                 hprof2transverse->Fill( ptleadjet, ptsumevent);
 867                 hpttransverse->Fill( ptpart );
 868               }
 869             else
 870               if ( ( dphipartljet[part] >= ( 2.*kPI/3.)) && ( dphipartljet[part] < kPI ) )
 871                 {
 872                   hprof1away->Fill( ptleadjet, fNPart );
 873                   hprof2away->Fill( ptleadjet, ptsumevent);
 874                   hptaway->Fill( ptpart );
 875                 }
 876         }
 877     }
 878
 879   delete [] dphipartljet;
 880   delete [] zpartljet;
 881   delete [] idxpartLJ;
 882
 883 }
 884
 885 //______________________________________________________________________________
 886 void AliCdfJetFinder::Clean()
 887 {
 888   // CLEANING SECTION
 889   for (  Int_t i = 0 ; i < fNPart ; i++ ){
 890     delete fVectParticle[i];
 891     fVectParticle[i] = 0;
 892   }
 893   delete [] fVectParticle;fVectParticle = 0;
 894
 895   for (  Int_t i = 0 ; i < fNJets ; i++ ){
 896     delete fVectJet[i];
 897     fVectJet[i] = 0;
 898   }
 899   delete [] fVectJet;fVectJet = 0;
 900
 901   delete [] fPtArray;fPtArray = 0;
 902   delete [] fIdxArray;fIdxArray = 0;
 903
 904   Reset();
 905
 906 }
 907