EMCAL/AliEMCALJetFinder.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 $Log$
  18 Revision 1.36  2003/01/15 19:05:44  morsch
  19 Updated selection in ReadFromTracks()
  20
  21 Revision 1.35  2003/01/15 04:59:38  morsch
  22 - TPC eff. from AliEMCALFast
  23 - Correction in PropagatePhi()
  24
  25 Revision 1.34  2003/01/14 10:50:18  alibrary
  26 Cleanup of STEER coding conventions
  27
  28 Revision 1.33  2003/01/10 10:48:19  morsch
  29 SetSamplingFraction() removed from constructor.
  30
  31 Revision 1.32  2003/01/10 10:26:40  morsch
  32 Sampling fraction initialized from geometry class.
  33
  34 Revision 1.31  2003/01/08 17:13:41  schutz
  35 added the HCAL section
  36
  37 Revision 1.30  2002/12/09 16:26:28  morsch
  38 - Nummber of particles per jet increased to 1000
  39 - Warning removed.
  40
  41 Revision 1.29  2002/11/21 17:01:40  alibrary
  42 Removing AliMCProcess and AliMC
  43
  44 Revision 1.28  2002/11/20 14:13:16  morsch
  45 - FindChargedJets() added.
  46 - Destructor corrected.
  47 - Geometry cuts taken from AliEMCALGeometry.
  48
  49 Revision 1.27  2002/11/15 17:39:10  morsch
  50 GetPythiaParticleName removed.
  51
  52 Revision 1.26  2002/10/14 14:55:35  hristov
  53 Merging the VirtualMC branch to the main development branch (HEAD)
  54
  55 Revision 1.20.4.3  2002/10/10 15:07:49  hristov
  56 Updating VirtualMC to v3-09-02
  57
  58 Revision 1.25  2002/09/13 10:24:57  morsch
  59 idem
  60
  61 Revision 1.24  2002/09/13 10:21:13  morsch
  62 No cast to AliMagFCM.
  63
  64 Revision 1.23  2002/06/27 09:24:26  morsch
  65 Uncomment the TH1::AddDirectory statement.
  66
  67 Revision 1.22  2002/05/22 13:48:43  morsch
  68 Pdg code added to track list.
  69
  70 Revision 1.21  2002/04/27 07:43:08  morsch
  71 Calculation of fDphi corrected (Renan Cabrera)
  72
  73 Revision 1.20  2002/03/12 01:06:23  pavlinov
  74 Testin output from generator
  75
  76 Revision 1.19  2002/02/27 00:46:33  pavlinov
  77 Added method FillFromParticles()
  78
  79 Revision 1.18  2002/02/21 08:48:59  morsch
  80 Correction in FillFromHitFlaggedTrack. (Jennifer Klay)
  81
  82 Revision 1.17  2002/02/14 08:52:53  morsch
  83 Major updates by Aleksei Pavlinov:
  84 FillFromPartons, FillFromTracks, jetfinder configuration.
  85
  86 Revision 1.16  2002/02/08 11:43:05  morsch
  87 SetOutputFileName(..) allows to specify an output file to which the
  88 reconstructed jets are written. If not set output goes to input file.
  89 Attention call Init() before processing.
  90
  91 Revision 1.15  2002/02/02 08:37:18  morsch
  92 Formula for rB corrected.
  93
  94 Revision 1.14  2002/02/01 08:55:30  morsch
  95 Fill map with Et and pT.
  96
  97 Revision 1.13  2002/01/31 09:37:36  morsch
  98 Geometry parameters in constructor and call of SetCellSize()
  99
 100 Revision 1.12  2002/01/23 13:40:23  morsch
 101 Fastidious debug print statement removed.
 102
 103 Revision 1.11  2002/01/22 17:25:47  morsch
 104 Some corrections for event mixing and bg event handling.
 105
 106 Revision 1.10  2002/01/22 10:31:44  morsch
 107 Some correction for bg mixing.
 108
 109 Revision 1.9  2002/01/21 16:28:42  morsch
 110 Correct sign of dphi.
 111
 112 Revision 1.8  2002/01/21 16:05:31  morsch
 113 - different phi-bin for hadron correction
 114 - provisions for background mixing.
 115
 116 Revision 1.7  2002/01/21 15:47:18  morsch
 117 Bending radius correctly in cm.
 118
 119 Revision 1.6  2002/01/21 12:53:50  morsch
 120 authors
 121
 122 Revision 1.5  2002/01/21 12:47:47  morsch
 123 Possibility to include K0long and neutrons.
 124
 125 Revision 1.4  2002/01/21 11:03:21  morsch
 126 Phi propagation introduced in FillFromTracks.
 127
 128 Revision 1.3  2002/01/18 05:07:56  morsch
 129 - hadronic correction
 130 - filling of digits
 131 - track selection upon EMCAL information
 132
 133 */
 134
 135 //*-- Authors: Andreas Morsch   (CERN)
 136 //*            J.L. Klay        (LBL)
 137 //*            Aleksei Pavlinov (WSU)
 138
 139 #include <stdio.h>
 140 // From root ...
 141 #include <TArrayF.h>
 142 #include <TAxis.h>
 143 #include <TBranchElement.h>
 144 #include <TCanvas.h>
 145 #include <TClonesArray.h>
 146 #include <TFile.h>
 147 #include <TH1.h>
 148 #include <TH2.h>
 149 #include <TList.h>
 150 #include <TPDGCode.h>
 151 #include <TPad.h>
 152 #include <TParticle.h>
 153 #include <TParticlePDG.h>
 154 #include <TPaveText.h>
 155 #include <TPythia6Calls.h>
 156 #include <TROOT.h>
 157 #include <TStyle.h>
 158 #include <TTree.h>
 159
 160 // From AliRoot ...
 161 #include "AliEMCAL.h"
 162 #include "AliEMCALDigit.h"
 163 #include "AliEMCALDigitizer.h"
 164 #include "AliEMCALFast.h"
 165 #include "AliEMCALGeometry.h"
 166 #include "AliEMCALHadronCorrection.h"
 167 #include "AliEMCALHit.h"
 168 #include "AliEMCALJetFinder.h"
 169 #include "AliEMCALJetMicroDst.h"
 170 #include "AliHeader.h"
 171 #include "AliMagF.h"
 172 #include "AliMagFCM.h"
 173 #include "AliRun.h"
 174
 175 // Interface to FORTRAN
 176 #include "Ecommon.h"
 177
 178
 179 ClassImp(AliEMCALJetFinder)
 180
 181 //____________________________________________________________________________
 182 AliEMCALJetFinder::AliEMCALJetFinder()
 183 {
 184 // Default constructor
 185     fJets             = 0;
 186     fNjets            = 0;
 187     fLego             = 0;
 188     fLegoB            = 0;
 189
 190     fTrackList        = 0;
 191     fPtT              = 0;
 192     fEtaT             = 0;
 193     fPhiT             = 0;
 194     fPdgT             = 0;
 195
 196     fTrackListB       = 0;
 197     fPtB              = 0;
 198     fEtaB             = 0;
 199     fPhiB             = 0;
 200     fPdgB             = 0;
 201
 202     fHCorrection      = 0;
 203     fHadronCorrector  = 0;
 204
 205     fWrite            = 0;
 206     fOutFileName      = 0;
 207     fOutFile          = 0;
 208     fInFile           = 0;
 209     fEvent            = 0;
 210
 211     SetParametersForBgSubtraction();
 212 }
 213
 214 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
 215     : TTask(name, title)
 216 {
 217 // Constructor
 218 // Title is used in method GetFileNameForParameters();
 219 //
 220     fJets  = new TClonesArray("AliEMCALJet",10000);
 221     fNjets = 0;
 222     for (Int_t i = 0; i < 30000; i++)
 223     {
 224         fEtCell[i]  = 0.;
 225         fEtaCell[i] = 0.;
 226         fPhiCell[i] = 0.;
 227     }
 228     fLego       = 0;
 229     fLegoB      = 0;
 230
 231     fTrackList  = 0;
 232     fPtT        = 0;
 233     fEtaT       = 0;
 234     fPhiT       = 0;
 235     fPdgT       = 0;
 236
 237     fTrackListB       = 0;
 238     fPtB        = 0;
 239     fEtaB       = 0;
 240     fPhiB       = 0;
 241     fPdgB       = 0;
 242
 243     fHCorrection      = 0;
 244     fHadronCorrector  = 0;
 245     fBackground       = 0;
 246     fWrite            = 0;
 247     fOutFileName      = 0;
 248     fOutFile          = 0;
 249     fInFile           = 0;
 250     fEvent            = 0;
 251 //
 252     SetPtCut();
 253     SetMomentumSmearing();
 254     SetEfficiencySim();
 255     SetDebug();
 256     SetHadronCorrection();
 257     SetIncludeK0andN();
 258
 259     SetParametersForBgSubtraction();
 260 }
 261
 262 void AliEMCALJetFinder::SetParametersForBgSubtraction
 263 (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
 264 {
 265 // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
 266 // at WSU Linux cluster - 11-feb-2002
 267 // These parameters must be tuned more carefull !!!
 268   SetMode(mode);
 269   SetMinMove(minMove);
 270   SetMaxMove(maxMove);
 271   SetPrecBg(precBg);
 272 }
 273
 274 //____________________________________________________________________________
 275 AliEMCALJetFinder::~AliEMCALJetFinder()
 276 {
 277 // Destructor
 278 //
 279     if (fJets){
 280         fJets->Delete();
 281         delete fJets;
 282     }
 283     delete fLego;
 284     delete fLegoB;
 285     delete fhLegoTracks;
 286     delete fhLegoEMCAL;
 287     delete fhLegoHadrCorr;
 288     delete fhEff;
 289     delete fhCellEt;
 290     delete fhCellEMCALEt;
 291     delete fhTrackPt;
 292     delete fhTrackPtBcut;
 293     delete fhChPartMultInTpc;
 294
 295     delete[] fTrackList;
 296     delete[] fPtT;
 297     delete[] fEtaT;
 298     delete[] fPhiT;
 299     delete[] fPdgT;
 300
 301     delete[] fTrackListB;
 302     delete[] fPtB;
 303     delete[] fEtaB;
 304     delete[] fPhiB;
 305     delete[] fPdgB;
 306 }
 307
 308 #ifndef WIN32
 309 # define jet_finder_ua1 jet_finder_ua1_
 310 # define hf1 hf1_
 311 # define type_of_call
 312
 313 #else
 314 # define jet_finder_ua1 JET_FINDER_UA1
 315 # define hf1 HF1
 316 # define type_of_call _stdcall
 317 #endif
 318
 319 extern "C" void type_of_call
 320 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
 321                Float_t  etc[30000],  Float_t etac[30000],
 322                Float_t  phic[30000],
 323                Float_t& min_move, Float_t& max_move, Int_t& mode,
 324                Float_t& prec_bg,  Int_t& ierror);
 325
 326 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
 327
 328
 329 void AliEMCALJetFinder::Init()
 330 {
 331 //
 332 // Geometry and I/O initialization
 333 //
 334 //
 335 //
 336 //  Get geometry parameters from EMCAL
 337 //
 338 //
 339 //  Geometry
 340     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 341     AliEMCALGeometry* geom =
 342     AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 343
 344 //    SetSamplingFraction(geom->GetSampling());
 345
 346     fNbinEta = geom->GetNZ();
 347     fNbinPhi = geom->GetNPhi();
 348     fPhiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
 349     fPhiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
 350     fEtaMin  = geom->GetArm1EtaMin();
 351     fEtaMax  = geom->GetArm1EtaMax();
 352     fDphi    = (fPhiMax-fPhiMin)/fNbinPhi;
 353     fDeta    = (fEtaMax-fEtaMin)/fNbinEta;
 354     fNtot    = fNbinPhi*fNbinEta;
 355     fWeightingMethod = kFALSE;
 356
 357 //
 358     SetCellSize(fDeta, fDphi);
 359 //
 360 //  I/O
 361     if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
 362 }
 363
 364 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
 365                              Float_t etac[30000], Float_t phic[30000],
 366                              Float_t min_move, Float_t max_move, Int_t mode,
 367                              Float_t prec_bg,  Int_t   ierror)
 368 {
 369 // Wrapper for fortran coded jet finder
 370 // Full argument list
 371     jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
 372                    min_move, max_move, mode, prec_bg, ierror);
 373     // Write result to output
 374     if(fWrite) WriteJets();
 375     fEvent++;
 376 }
 377
 378 void AliEMCALJetFinder::Find()
 379 {
 380 // Wrapper for fortran coded jet finder using member data for
 381 // argument list
 382
 383     Float_t min_move = fMinMove;
 384     Float_t max_move = fMaxMove;
 385     Int_t   mode     = fMode;
 386     Float_t prec_bg  = fPrecBg;
 387     Int_t   ierror;
 388     ResetJets(); // 4-feb-2002 by PAI
 389
 390     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
 391                    min_move, max_move, mode, prec_bg, ierror);
 392     fError = ierror;
 393     // Write result to output
 394     Int_t njet = Njets();
 395
 396     for (Int_t nj=0; nj<njet; nj++)
 397     {
 398         if (fWeightingMethod)
 399         {
 400           fJetT[nj] = new AliEMCALJet(WeightedJetEnergy( JetEtaW(nj),JetPhiW(nj) ),
 401                                       JetPhiW(nj),
 402                                       JetEtaW(nj));
 403
 404         }else{
 405
 406         fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
 407                                     JetPhiW(nj),
 408                                     JetEtaW(nj));
 409         }
 410         fJetT[nj]->SetIsWeightedEnergy(fWeightingMethod);
 411         fJetT[nj]->SetEMCALEnergy( EMCALConeEnergy(JetEtaW(nj),JetPhiW(nj)) );
 412         fJetT[nj]->SetTrackEnergy(TrackConeEnergy( JetEtaW(nj),JetPhiW(nj)  ));
 413         fJetT[nj]->SetHCEnergy(HCConeEnergy( JetEtaW(nj),JetPhiW(nj)  ));
 414
 415     }
 416
 417     FindTracksInJetCone();
 418     if(fWrite) WriteJets();
 419     fEvent++;
 420 }
 421
 422
 423 Float_t AliEMCALJetFinder::EMCALConeEnergy(Float_t eta, Float_t phi)
 424 {
 425 Float_t newenergy = 0.0;
 426 Float_t bineta,binphi;
 427 TAxis *x = fhLegoEMCAL->GetXaxis();
 428 TAxis *y = fhLegoEMCAL->GetYaxis();
 429 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 430         {
 431    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 432       {
 433          bineta = x->GetBinCenter(i);
 434          binphi = y->GetBinCenter(j);
 435        if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 436       {
 437            newenergy += fhLegoEMCAL->GetBinContent(i,j);
 438       }
 439     }
 440 }
 441 return newenergy;
 442 }
 443
 444 Float_t AliEMCALJetFinder::TrackConeEnergy(Float_t eta, Float_t phi)
 445 {
 446 Float_t newenergy = 0.0;
 447 Float_t bineta,binphi;
 448 TAxis *x = fhLegoTracks->GetXaxis();
 449 TAxis *y = fhLegoTracks->GetYaxis();
 450 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 451    {
 452    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 453    {
 454     bineta = x->GetBinCenter(i);
 455     binphi = y->GetBinCenter(j);
 456     if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 457       {
 458           newenergy += fhLegoTracks->GetBinContent(i,j);
 459        }
 460     }
 461 }
 462 return newenergy;
 463 }
 464
 465 Float_t AliEMCALJetFinder::HCConeEnergy(Float_t eta, Float_t phi)
 466 {
 467 //Float_t newenergy = 0.0;
 468 //Float_t bineta,binphi;
 469 //TAxis *x = fhLegoTracks->GetXaxis();
 470 //TAxis *y = fhLegoTracks->GetYaxis();
 471 //for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 472 //   {
 473 //      for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 474 //      {
 475 //        bineta = x->GetBinCenter(i);
 476 //        binphi = y->GetBinCenter(j);
 477 //        if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 478 //        {
 479 //            newenergy += fhLegoTracks->GetBinContent(i,j);
 480 //         }
 481 //    }
 482 //}
 483 //return newenergy;
 484
 485 return 0.0;
 486
 487 }
 488
 489
 490
 491 Float_t AliEMCALJetFinder::WeightedJetEnergy(Float_t eta, Float_t phi)
 492 {
 493
 494
 495 Float_t newenergy = 0.0;
 496 Float_t bineta,binphi;
 497 TAxis *x = fhLegoEMCAL->GetXaxis();
 498 TAxis *y = fhLegoEMCAL->GetYaxis();
 499
 500
 501 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 502 {
 503    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 504    {
 505       bineta = x->GetBinCenter(i);
 506       binphi = y->GetBinCenter(j);
 507       if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 508       {
 509           newenergy += (fEMCALWeight)* fhLegoEMCAL->GetBinContent(i,j) + (fTrackWeight)* fhLegoTracks->GetBinContent(i,j);
 510       }
 511     }
 512 }
 513
 514 return newenergy;
 515
 516 }
 517
 518
 519 Int_t AliEMCALJetFinder::Njets()
 520 {
 521 // Get number of reconstructed jets
 522     return EMCALJETS.njet;
 523 }
 524
 525 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
 526 {
 527 // Get reconstructed jet energy
 528     return EMCALJETS.etj[i];
 529 }
 530
 531 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
 532 {
 533 // Get reconstructed jet phi from leading particle
 534     return EMCALJETS.phij[0][i];
 535 }
 536
 537 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
 538 {
 539 // Get reconstructed jet phi from weighting
 540     return EMCALJETS.phij[1][i];
 541 }
 542
 543 Float_t  AliEMCALJetFinder::JetEtaL(Int_t i)
 544 {
 545 // Get reconstructed jet eta from leading particles
 546     return EMCALJETS.etaj[0][i];
 547 }
 548
 549
 550 Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)
 551 {
 552 // Get reconstructed jet eta from weighting
 553     return EMCALJETS.etaj[1][i];
 554 }
 555
 556 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
 557 {
 558 // Set grid cell size
 559     EMCALCELLGEO.etaCellSize = eta;
 560     EMCALCELLGEO.phiCellSize = phi;
 561 }
 562
 563 void AliEMCALJetFinder::SetConeRadius(Float_t par)
 564 {
 565 // Set jet cone radius
 566     EMCALJETPARAM.coneRad = par;
 567     fConeRadius = par;
 568 }
 569
 570 void AliEMCALJetFinder::SetEtSeed(Float_t par)
 571 {
 572 // Set et cut for seeds
 573     EMCALJETPARAM.etSeed = par;
 574     fEtSeed = par;
 575 }
 576
 577 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
 578 {
 579 // Set minimum jet et
 580     EMCALJETPARAM.ejMin = par;
 581     fMinJetEt = par;
 582 }
 583
 584 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
 585 {
 586 // Set et cut per cell
 587     EMCALJETPARAM.etMin = par;
 588     fMinCellEt = par;
 589 }
 590
 591 void AliEMCALJetFinder::SetPtCut(Float_t par)
 592 {
 593 // Set pT cut on charged tracks
 594     fPtCut = par;
 595 }
 596
 597
 598 void AliEMCALJetFinder::Test()
 599 {
 600 //
 601 // Test the finder call
 602 //
 603     const Int_t nmax = 30000;
 604     Int_t ncell      = 10;
 605     Int_t ncell_tot  = 100;
 606
 607     Float_t etc[nmax];
 608     Float_t etac[nmax];
 609     Float_t phic[nmax];
 610     Float_t min_move = 0;
 611     Float_t max_move = 0;
 612     Int_t   mode = 0;
 613     Float_t prec_bg = 0;
 614     Int_t   ierror = 0;
 615
 616
 617     Find(ncell, ncell_tot, etc, etac, phic,
 618          min_move, max_move, mode, prec_bg, ierror);
 619
 620 }
 621
 622 //
 623 //  I/O
 624 //
 625
 626 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
 627 {
 628     //
 629     // Add a jet
 630     //
 631     TClonesArray &lrawcl = *fJets;
 632     new(lrawcl[fNjets++]) AliEMCALJet(jet);
 633 }
 634
 635 void AliEMCALJetFinder::ResetJets()
 636 {
 637     //
 638     // Reset Jet List
 639     //
 640     fJets->Clear();
 641     fNjets = 0;
 642 }
 643
 644 void AliEMCALJetFinder::WriteJets()
 645 {
 646 //
 647 // Add all jets to the list
 648 //
 649     const Int_t kBufferSize = 4000;
 650     const char* file = 0;
 651
 652     Int_t njet = Njets();
 653
 654     for (Int_t nj = 0; nj < njet; nj++)
 655     {
 656         AddJet(*fJetT[nj]);
 657         delete fJetT[nj];
 658     }
 659
 660 // I/O
 661     if (!fOutFileName) {
 662 //
 663 // output written to input file
 664 //
 665         AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
 666         TTree* pK = gAlice->TreeK();
 667         file = (pK->GetCurrentFile())->GetName();
 668         if (fDebug > 1)
 669             printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
 670         if (fJets && gAlice->TreeR()) {
 671             pEMCAL->MakeBranchInTree(gAlice->TreeR(),
 672                                      "EMCALJets",
 673                                      &fJets,
 674                                      kBufferSize,
 675                                      file);
 676         }
 677         Int_t nev = gAlice->GetHeader()->GetEvent();
 678         gAlice->TreeR()->Fill();
 679         char hname[30];
 680         sprintf(hname,"TreeR%d", nev);
 681         gAlice->TreeR()->Write(hname);
 682         gAlice->TreeR()->Reset();
 683     } else {
 684 //
 685 // Output written to user specified output file
 686 //
 687         TTree* pK = gAlice->TreeK();
 688         fInFile  = pK->GetCurrentFile();
 689
 690         fOutFile->cd();
 691         char hname[30];
 692         sprintf(hname,"TreeR%d", fEvent);
 693         TTree* treeJ = new TTree(hname, "EMCALJets");
 694         treeJ->Branch("EMCALJets", &fJets, kBufferSize);
 695         treeJ->Fill();
 696         treeJ->Write(hname);
 697         fInFile->cd();
 698     }
 699     ResetJets();
 700 }
 701
 702 void AliEMCALJetFinder::BookLego()
 703 {
 704 //
 705 //  Book histo for discretisation
 706 //
 707
 708 //
 709 //  Don't add histos to the current directory
 710     if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
 711
 712     TH2::AddDirectory(0);
 713     TH1::AddDirectory(0);
 714     gROOT->cd();
 715 //
 716 //  Signal map
 717     fLego = new TH2F("legoH","eta-phi",
 718                            fNbinEta, fEtaMin, fEtaMax,
 719                            fNbinPhi, fPhiMin, fPhiMax);
 720 //
 721 //  Background map
 722     fLegoB = new TH2F("legoB","eta-phi for BG event",
 723                            fNbinEta, fEtaMin, fEtaMax,
 724                            fNbinPhi, fPhiMin, fPhiMax);
 725
 726 //  Tracks map
 727     fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks",
 728     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 729 //  EMCAL map
 730     fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL",
 731     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 732 //  Hadron correction map
 733     fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction",
 734     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 735 //  Hists. for tuning jet finder
 736     fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.);
 737
 738     TArrayF eTmp(1101);
 739     eTmp[0] = 0.0;
 740     for(Int_t i=1; i<=1000; i++)      eTmp[i] = 0.1*i; // step 100 mev
 741     for(Int_t i=1001; i<=1100; i++)   eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV
 742
 743     fhCellEt  = new TH1F("hCellEt","Cell E_{T} from fLego",
 744     eTmp.GetSize()-1, eTmp.GetArray());
 745     fhCellEMCALEt  = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself",
 746     eTmp.GetSize()-1, eTmp.GetArray());
 747     fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ",
 748     eTmp.GetSize()-1, eTmp.GetArray());
 749     fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
 750     eTmp.GetSize()-1, eTmp.GetArray());
 751
 752     fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
 753     "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
 754
 755             //! first canvas for drawing
 756     fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kTRUE);
 757 }
 758
 759 void AliEMCALJetFinder::DumpLego()
 760 {
 761 //
 762 // Dump lego histo into array
 763 //
 764     fNcell = 0;
 765     TAxis* Xaxis = fLego->GetXaxis();
 766     TAxis* Yaxis = fLego->GetYaxis();
 767     //    fhCellEt->Clear();
 768     Float_t e, eH;
 769     for (Int_t i = 1; i <= fNbinEta; i++) {
 770         for (Int_t j = 1; j <= fNbinPhi; j++) {
 771             e = fLego->GetBinContent(i,j);
 772             if (e > 0.0) {
 773               Float_t eta  = Xaxis->GetBinCenter(i);
 774               Float_t phi  = Yaxis->GetBinCenter(j);
 775               fEtCell[fNcell]  = e;
 776               fEtaCell[fNcell] = eta;
 777               fPhiCell[fNcell] = phi;
 778               fNcell++;
 779               fhCellEt->Fill(e);
 780             }
 781             if(fhLegoEMCAL) {
 782               eH = fhLegoEMCAL->GetBinContent(i,j);
 783               if(eH > 0.0) fhCellEMCALEt->Fill(eH);
 784             }
 785         } // phi
 786     } // eta
 787 //  today
 788 //    fNcell--;
 789 }
 790
 791 void AliEMCALJetFinder::ResetMap()
 792 {
 793 //
 794 // Reset eta-phi array
 795
 796     for (Int_t i=0; i<30000; i++)
 797     {
 798         fEtCell[i]  = 0.;
 799         fEtaCell[i] = 0.;
 800         fPhiCell[i] = 0.;
 801     }
 802 }
 803
 804
 805 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 806 {
 807 //
 808 // Fill Cells with track information
 809 //
 810     if (fDebug >= 2)
 811     printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
 812     fNChTpc = 0;
 813
 814     ResetMap();
 815
 816     if (!fLego) BookLego();
 817 // Reset
 818     if (flag == 0) fLego->Reset();
 819 //
 820 // Access particle information
 821     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 822     Int_t ntr   = (gAlice->GetHeader())->GetNtrack();
 823     printf(" : #primary particles %i # tracks %i \n", npart, ntr);
 824
 825 // Create track list
 826 //
 827 // 0: not selected
 828 // 1: selected for jet finding
 829 // 2: inside jet
 830 // ....
 831     if (fTrackList) delete[] fTrackList;
 832     if (fPtT)       delete[] fPtT;
 833     if (fEtaT)      delete[] fEtaT;
 834     if (fPhiT)      delete[] fPhiT;
 835     if (fPdgT)      delete[] fPdgT;
 836
 837     fTrackList = new Int_t  [npart];
 838     fPtT       = new Float_t[npart];
 839     fEtaT      = new Float_t[npart];
 840     fPhiT      = new Float_t[npart];
 841     fPdgT      = new Int_t[npart];
 842
 843     fNt   = npart;
 844     fNtS  = 0;
 845     Float_t chTmp=0.0; // charge of current particle
 846     //    Int_t idGeant;
 847
 848     // this is for Pythia ??
 849     for (Int_t part = 0; part < npart; part++) {
 850         TParticle *MPart = gAlice->Particle(part);
 851         Int_t mpart   = MPart->GetPdgCode();
 852         Int_t child1  = MPart->GetFirstDaughter();
 853         Float_t pT    = MPart->Pt();
 854         Float_t p     = MPart->P();
 855         Float_t phi   = MPart->Phi();
 856         Float_t eta   = -100.;
 857         if(pT > 0.001) eta   = MPart->Eta();
 858         Float_t theta = MPart->Theta();
 859         if  (fDebug>=2) {
 860            printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
 861            part, mpart, child1, MPart->GetLastDaughter(), MPart->GetStatusCode());
 862         }
 863
 864         if (fDebug >= 2) {
 865             if (part == 6 || part == 7)
 866             {
 867                 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
 868                        part-5, pT, eta, phi);
 869             }
 870         }
 871
 872         fTrackList[part] = 0;
 873         fPtT[part]       = pT; // must be change after correction for resolution !!!
 874         fEtaT[part]      = eta;
 875         fPhiT[part]      = phi;
 876         fPdgT[part]      = mpart;
 877
 878 //      if (part < 2) continue;
 879         if (MPart->GetStatusCode() != 1) continue;
 880
 881         TParticlePDG* pdgP = 0;
 882 // charged or neutral
 883         pdgP  = MPart->GetPDG();
 884         chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
 885
 886         if (ich == 0) {
 887             if (chTmp == 0) {
 888                 if (!fK0N) {
 889                     continue;
 890                 } else {
 891                     if (mpart != kNeutron    &&
 892                         mpart != kNeutronBar &&
 893                         mpart != kK0Long) continue;
 894                 }
 895             }
 896         } else if (ich == 2) {
 897           if (mpart == kNeutron    ||
 898               mpart == kNeutronBar ||
 899               mpart == kK0Long) continue;
 900         }
 901
 902         if (TMath::Abs(eta)<=0.9) fNChTpc++;
 903 // final state only
 904         if (child1 >= 0 && child1 < npart) continue;
 905 // acceptance cut
 906         if (eta > fEtaMax || eta < fEtaMin)    continue;
 907         if (phi > fPhiMax || phi < fPhiMin )   continue;
 908 //
 909         if (fDebug >= 3)
 910         printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
 911         part, mpart, child1, eta, phi, pT, chTmp);
 912 //
 913 //
 914 // Momentum smearing goes here ...
 915 //
 916         fhTrackPt->Fill(pT);
 917         Float_t pw;
 918         if (fSmear && TMath::Abs(chTmp)) {
 919             pw = AliEMCALFast::SmearMomentum(1,p);
 920             // p changed - take into account when calculate pT,
 921             // pz and so on ..  ?
 922             pT = (pw/p) * pT;
 923             if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
 924             p  = pw;
 925         }
 926 //
 927 // Tracking Efficiency and TPC acceptance goes here ...
 928         Float_t eff;
 929         if (fEffic && TMath::Abs(chTmp)) {
 930             eff = 0.9; //  AliEMCALFast::Efficiency(2,p);
 931             if(fhEff) fhEff->Fill(p, eff);
 932             if (AliEMCALFast::RandomReject(eff)) {
 933                 if(fDebug >= 5) printf(" reject due to unefficiency ");
 934                 continue;
 935             }
 936         }
 937 //
 938 // Correction of Hadronic Energy goes here ...
 939 //
 940 //
 941 // phi propagation for hadronic correction
 942
 943         Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
 944         Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
 945         if(TMath::Abs(chTmp)) {
 946             // hadr. correction only for charge particle
 947             dphi  = PropagatePhi(pT, chTmp, curls);
 948             phiHC = phi + dphi;
 949             if (fDebug >= 6) {
 950                 printf("\n Delta phi %f pT %f ", dphi, pT);
 951                 if (curls) printf("\n !! Track is curling");
 952             }
 953             if(!curls) fhTrackPtBcut->Fill(pT);
 954
 955             if (fHCorrection && !curls) {
 956                 if (!fHadronCorrector)
 957                     Fatal("AliEMCALJetFinder",
 958                           "Hadronic energy correction required but not defined !");
 959
 960                 dpH    = fHadronCorrector->GetEnergy(p, eta, 7);
 961                 eTdpH  = dpH*TMath::Sin(theta);
 962
 963                 if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n",
 964                                         phi, phiHC, -eTdpH); // correction is negative
 965                 Int_t xbin,ybin;
 966                 xbin = fLego->GetXaxis()->FindBin(eta);
 967                 ybin = fLego->GetYaxis()->FindBin(phiHC);
 968                 cout <<"Hadron Correction affected bin - contents before correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
 969                 fLego->Fill(eta, phi, -fSamplingF*eTdpH );
 970                 cout <<"Hadron Correction affected bin - contents after  correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
 971                 fhLegoHadrCorr->Fill(eta, phi, fSamplingF*eTdpH);
 972             }
 973         }
 974 //
 975 //  More to do ? Just think about it !
 976 //
 977         if (phi > fPhiMax || phi < fPhiMin )   continue;
 978
 979         if(TMath::Abs(chTmp) ) { // charge particle
 980             if (pT > fPtCut && !curls) {
 981                 if (fDebug >= 8) printf("Charge :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
 982                                         eta , phi, pT, fNtS);
 983                 fLego->Fill(eta, phi, pT);
 984                 fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
 985                 fTrackList[part] = 1;
 986                 fNtS++;
 987             }
 988         } else if(ich > 0 || fK0N) {
 989             // case of n, nbar and K0L
 990             if (fDebug >= 9) printf("Neutral :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
 991                                     eta , phi, pT, fNtS);
 992             fLego->Fill(eta, phi, pT);
 993             fTrackList[part] = 1;
 994             fNtS++;
 995         }
 996     } // primary loop
 997     Float_t etsum = 0.;
 998     for(Int_t i=0; i<fLego->GetSize(); i++) {
 999         Float_t etc =  (*fLego)[i];
1000         if (etc > fMinCellEt) etsum += etc;
1001     }
1002
1003     printf("\nFillFromTracks: Sum above threshold %f %f \n \n", fMinCellEt, etsum);
1004
1005     DumpLego();
1006 }
1007
1008 void AliEMCALJetFinder::FillFromHits(Int_t flag)
1009 {
1010 //
1011 // Fill Cells with hit information
1012 //
1013 //
1014     if (fDebug >= 2)
1015     printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
1016
1017     ResetMap();
1018
1019     if (!fLego) BookLego();
1020 //  Reset eta-phi maps if needed
1021     if (flag == 0) { // default behavior
1022       fLego->Reset();
1023       fhLegoTracks->Reset();
1024       fhLegoEMCAL->Reset();
1025       fhLegoHadrCorr->Reset();
1026     }
1027 //  Initialize from background event if available
1028 //
1029 // Access hit information
1030     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1031     TTree *treeH = gAlice->TreeH();
1032     Int_t ntracks = (Int_t) treeH->GetEntries();
1033 //
1034 //   Loop over tracks
1035 //
1036     Int_t nbytes = 0;
1037     Double_t etH = 0.0;
1038
1039     for (Int_t track=0; track<ntracks;track++) {
1040         gAlice->ResetHits();
1041         nbytes += treeH->GetEvent(track);
1042 //
1043 //  Loop over hits
1044 //
1045         for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1046             mHit;
1047             mHit=(AliEMCALHit*) pEMCAL->NextHit())
1048         {
1049             Float_t x      =    mHit->X();         // x-pos of hit
1050             Float_t y      =    mHit->Y();         // y-pos
1051             Float_t z      =    mHit->Z();         // z-pos
1052             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
1053
1054             Float_t r      =    TMath::Sqrt(x*x+y*y);
1055             Float_t theta  =    TMath::ATan2(r,z);
1056             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
1057             Float_t phi    =    TMath::ATan2(y,x);
1058
1059             if (fDebug >= 11) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
1060 //          printf("\n Hit %f %f %f %f", x, y, z, eloss);
1061
1062             etH = fSamplingF*eloss*TMath::Sin(theta);
1063             fLego->Fill(eta, phi, etH);
1064             //      fhLegoEMCAL->Fill(eta, phi, etH);
1065         } // Hit Loop
1066     } // Track Loop
1067     // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
1068     Float_t etsum = 0;
1069
1070     for(Int_t i=0; i<fLego->GetSize(); i++) {
1071         (*fhLegoEMCAL)[i] = (*fLego)[i];
1072         Float_t etc =  (*fLego)[i];
1073         if (etc > fMinCellEt) etsum += etc;
1074     }
1075
1076     printf("\nFillFromHits: Sum above threshold %f %f \n \n", fMinCellEt, etsum);
1077
1078     //    DumpLego(); ??
1079
1080 }
1081
1082 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
1083 {
1084 //
1085 // Fill Cells with digit information
1086 //
1087 //
1088     ResetMap();
1089
1090     if (!fLego) BookLego();
1091     if (flag == 0) fLego->Reset();
1092     Int_t nbytes;
1093
1094
1095 //
1096 //  Connect digits
1097 //
1098     TClonesArray* digs      = new TClonesArray("AliEMCALDigit", 10000);
1099     TTree *treeD = gAlice->TreeD();
1100     TBranchElement* branchDg = (TBranchElement*)
1101         treeD->GetBranch("EMCAL");
1102
1103     if (!branchDg) Fatal("AliEMCALJetFinder",
1104                          "Reading digits requested but no digits in file !");
1105
1106     branchDg->SetAddress(&digs);
1107     Int_t nent = (Int_t) branchDg->GetEntries();
1108 //
1109 //  Connect digitizer
1110 //
1111     AliEMCALDigitizer* digr = new AliEMCALDigitizer();
1112     TBranchElement* branchDr = (TBranchElement*)
1113         treeD->GetBranch("AliEMCALDigitizer");
1114     branchDr->SetAddress(&digr);
1115 //
1116 //
1117     nbytes += branchDg->GetEntry(0);
1118     nbytes += branchDr->GetEntry(0);
1119 //
1120 //  Get digitizer parameters
1121     Float_t preADCped = digr->GetPREpedestal();
1122     Float_t preADCcha = digr->GetPREchannel();
1123     Float_t ecADCped  = digr->GetECpedestal();
1124     Float_t ecADCcha  = digr->GetECchannel();
1125     Float_t hcADCped  = digr->GetHCpedestal();
1126     Float_t hcADCcha  = digr->GetHCchannel();
1127
1128     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1129     AliEMCALGeometry* geom =
1130         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
1131
1132     if (fDebug) {
1133         Int_t   ndig = digs->GetEntries();
1134         Info("FillFromDigits","Number of Digits: %d %d\n Parameters: PRE : %f %f EC: %f %f HC: %f %f\n Geometry: %d %d",
1135              ndig, nent, preADCped, preADCcha, ecADCped, ecADCcha, hcADCped, hcADCcha, geom->GetNEta(), geom->GetNPhi());
1136     }
1137
1138 //
1139 //  Loop over digits
1140     AliEMCALDigit* sdg;
1141     TIter next(digs);
1142     while ((sdg = (AliEMCALDigit*)(next())))
1143     {
1144         Double_t pedestal = 0.;
1145         Double_t channel  = 0.;
1146         if (geom->IsInPRE(sdg->GetId())) {
1147           pedestal = preADCped;
1148           channel  = preADCcha;
1149         }
1150         else if (geom->IsInECAL(sdg->GetId()))  {
1151           pedestal = ecADCped;
1152           channel  = ecADCcha;
1153         }
1154         else if (geom->IsInHCAL(sdg->GetId()))  {
1155           pedestal = hcADCped;
1156           channel  = hcADCcha;
1157         }
1158         else
1159           Fatal("FillFromDigits", "unexpected digit is number!") ;
1160
1161         Float_t eta = sdg->GetEta();
1162         Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
1163         Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
1164
1165         if (fDebug > 1)
1166           Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
1167                            eta, phi, amp, sdg->GetAmp());
1168
1169         fLego->Fill(eta, phi, fSamplingF*amp);
1170     } // digit loop
1171 //
1172 //  Dump lego hist
1173     DumpLego();
1174 }
1175
1176
1177 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
1178 {
1179 //
1180 // Fill Cells with hit information
1181 //
1182 //
1183     ResetMap();
1184
1185     if (!fLego) BookLego();
1186 // Reset
1187     if (flag == 0) fLego->Reset();
1188
1189 // Flag charged tracks passing through TPC which
1190 // are associated to EMCAL Hits
1191     BuildTrackFlagTable();
1192
1193 //
1194 // Access particle information
1195     TTree *treeK = gAlice->TreeK();
1196     Int_t ntracks = (Int_t) treeK->GetEntries();
1197
1198     if (fPtT)       delete[] fPtT;
1199     if (fEtaT)      delete[] fEtaT;
1200     if (fPhiT)      delete[] fPhiT;
1201     if (fPdgT)      delete[] fPdgT;
1202
1203     fPtT       = new Float_t[ntracks];
1204     fEtaT      = new Float_t[ntracks];
1205     fPhiT      = new Float_t[ntracks];
1206     fPdgT      = new Int_t[ntracks];
1207
1208     fNt   = ntracks;
1209     fNtS  = 0;
1210
1211     for (Int_t track = 0; track < ntracks; track++) {
1212           TParticle *MPart = gAlice->Particle(track);
1213           Float_t pT    = MPart->Pt();
1214           Float_t phi   = MPart->Phi();
1215           Float_t eta   = MPart->Eta();
1216
1217         if(fTrackList[track]) {
1218           fPtT[track]       = pT;
1219           fEtaT[track]      = eta;
1220           fPhiT[track]      = phi;
1221           fPdgT[track]      = MPart->GetPdgCode();
1222
1223           if (track < 2) continue;      //Colliding particles?
1224           if (pT == 0 || pT < fPtCut) continue;
1225           fNtS++;
1226           fLego->Fill(eta, phi, pT);
1227         }
1228       } // track loop
1229     DumpLego();
1230 }
1231
1232 void AliEMCALJetFinder::FillFromParticles()
1233 {
1234 // 26-feb-2002 PAI - for checking all chain
1235 // Work on particles level; accept all particle (not neutrino )
1236
1237     Double_t PX=0, PY=0, PZ=0, E=0; // checking conservation law
1238     fNChTpc = 0;
1239
1240     ResetMap();
1241     if (!fLego) BookLego();
1242     fLego->Reset();
1243 //
1244 // Access particles information
1245     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1246     if (fDebug >= 2 || npart<=0) {
1247        printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
1248        if(npart<=0) return;
1249     }
1250     fNt   = npart;
1251     fNtS  = 0;
1252     RearrangeParticlesMemory(npart);
1253
1254 //  Go through the particles
1255     Int_t mpart, child1, child2, geantPdg;
1256     Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
1257     TParticle *MPart=0;
1258     for (Int_t part = 0; part < npart; part++) {
1259
1260         fTrackList[part] = 0;
1261
1262         MPart   = gAlice->Particle(part);
1263         mpart   = MPart->GetPdgCode();
1264         child1  = MPart->GetFirstDaughter();
1265         child2  = MPart->GetLastDaughter();
1266         pT      = MPart->Pt();
1267         phi     = MPart->Phi();
1268         eta     = MPart->Eta();
1269
1270         px      = MPart->Px();
1271         py      = MPart->Py();
1272         pz      = MPart->Pz();
1273         e       = MPart->Energy();
1274
1275 // see pyedit in Pythia's text
1276         geantPdg = mpart;
1277         if (IsThisPartonsOrDiQuark(mpart)) continue;
1278         printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n",
1279         part, mpart, geantPdg, px, py, pz, e, child1, child2);
1280
1281 //  exclude partons (21 - gluon, 92 - string)
1282
1283
1284 // exclude neutrinous also ??
1285         if (fDebug >= 11 && pT>0.01)
1286         printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
1287         part, mpart, eta, phi, pT);
1288
1289         fPtT[part]       = pT;
1290         fEtaT[part]      = eta;
1291         fPhiT[part]      = phi;
1292         fPdgT[part]      = mpart;
1293         fNtS++;
1294
1295 // final state only
1296         if (child1 >= 0 && child1 < npart) continue;
1297
1298         //        printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n",
1299         //        part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
1300         PX += px;
1301         PY += py;
1302         PZ += pz;
1303         E  += e;
1304
1305
1306         if (TMath::Abs(eta) <= 0.9) fNChTpc++;
1307 // acceptance cut
1308         if (eta > fEtaMax || eta < fEtaMin)    continue;
1309         if (phi > fPhiMax || phi < fPhiMin )   continue;
1310 //
1311         if(fK0N==0 ) { // exclude neutral hadrons
1312           if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue;
1313         }
1314         fTrackList[part] = 1;
1315         fLego->Fill(eta, phi, pT);
1316
1317     } // primary loop
1318     printf("\n                PX %8.2f  PY %8.2f  PZ %8.2f  E %8.2f \n",
1319     PX, PY, PZ, E);
1320     DumpLego();
1321     if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
1322 }
1323
1324 void AliEMCALJetFinder::FillFromPartons()
1325 {
1326 // function under construction - 13-feb-2002 PAI
1327
1328     if (fDebug >= 2)
1329     printf("\n AliEMCALJetFinder::FillFromPartons()\n");
1330     //
1331
1332     ResetMap();
1333     if (!fLego) BookLego();
1334     fLego->Reset();
1335 //
1336 // Access particle information
1337     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1338     if (fDebug >= 2 || npart<=0)
1339     printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
1340     fNt   = 0; // for FindTracksInJetCone
1341     fNtS  = 0;
1342
1343 //  Go through the partons
1344     Int_t statusCode=0;
1345     for (Int_t part = 8; part < npart; part++) {
1346         TParticle *MPart = gAlice->Particle(part);
1347         Int_t mpart   = MPart->GetPdgCode();
1348         //      Int_t child1  = MPart->GetFirstDaughter();
1349         Float_t pT    = MPart->Pt();
1350         //      Float_t p     = MPart->P();
1351         Float_t phi   = MPart->Phi();
1352         Float_t eta   = MPart->Eta();
1353         //      Float_t theta = MPart->Theta();
1354         statusCode    = MPart->GetStatusCode();
1355
1356 // accept partons (21 - gluon, 92 - string)
1357         if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
1358         if (fDebug > 1 && pT>0.01)
1359         printf("\n part:%5d mpart %5d status  %5d eta %8.2f phi %8.2f pT %8.2f ",
1360         part, mpart, statusCode, eta, phi, pT);
1361         //      if (fDebug >= 3) MPart->Print();
1362 // accept partons before fragmentation - p.57 in Pythia manual
1363 //        if(statusCode != 1) continue;
1364 // acceptance cut
1365         if (eta > fEtaMax || eta < fEtaMin)    continue;
1366         if (phi > fPhiMax || phi < fPhiMin )   continue;
1367 // final state only
1368 //      if (child1 >= 0 && child1 < npart) continue;
1369 //
1370 //
1371         fLego->Fill(eta, phi, pT);
1372
1373     } // primary loop
1374     DumpLego();
1375 }
1376
1377 void AliEMCALJetFinder::BuildTrackFlagTable() {
1378
1379 // Method to generate a lookup table for TreeK particles
1380 // which are linked to hits in the EMCAL
1381 //
1382 // --Author: J.L. Klay
1383 //
1384 // Access hit information
1385     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1386
1387     TTree *TK = gAlice->TreeK();                // Get the number of entries in the kine tree
1388     Int_t nKTrks = (Int_t) TK->GetEntries();    // (Number of particles created somewhere)
1389
1390     if(fTrackList) delete[] fTrackList;         //Make sure we get rid of the old one
1391     fTrackList = new Int_t[nKTrks];             //before generating a new one
1392
1393     for (Int_t i = 0; i < nKTrks; i++) {        //Initialize members to 0
1394         fTrackList[i] = 0;
1395     }
1396
1397     TTree *treeH = gAlice->TreeH();
1398     Int_t ntracks = (Int_t) treeH->GetEntries();
1399 //
1400 //   Loop over tracks
1401 //
1402     Int_t nbytes = 0;
1403
1404     for (Int_t track=0; track<ntracks;track++) {
1405         gAlice->ResetHits();
1406         nbytes += treeH->GetEvent(track);
1407         if (pEMCAL)  {
1408
1409 //
1410 //  Loop over hits
1411 //
1412             for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1413                 mHit;
1414                 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1415             {
1416                 Int_t   iTrk   =    mHit->Track();       // track number
1417                 Int_t   idprim =    mHit->GetPrimary();  // primary particle
1418
1419                 //Determine the origin point of this particle - it made a hit in the EMCAL
1420                 TParticle         *trkPart = gAlice->Particle(iTrk);
1421                 TParticlePDG  *trkPdg  = trkPart->GetPDG();
1422                 Int_t      trkCode = trkPart->GetPdgCode();
1423                 Double_t           trkChg;
1424                 if (trkCode < 10000) {          //Big Ions cause problems for
1425                     trkChg  = trkPdg->Charge(); //this function.  Since they aren't
1426                 } else {                                //likely to make it very far, set
1427                     trkChg  = 0.0;                      //their charge to 0 for the Flag test
1428                 }
1429                 Float_t            vxTrk   = trkPart->Vx();
1430                 Float_t            vyTrk   = trkPart->Vy();
1431                 Float_t            vrTrk   = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
1432                 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1433
1434                 //Loop through the ancestry of the EMCAL entrance particles
1435                 Int_t ancestor = trkPart->GetFirstMother();  //Get track's Mother
1436                 while (ancestor != -1) {
1437                     TParticle     *ancPart = gAlice->Particle(ancestor);  //get particle info on ancestor
1438                     TParticlePDG  *ancPdg  = ancPart->GetPDG();
1439                     Int_t             ancCode = ancPart->GetPdgCode();
1440                     Double_t       ancChg;
1441                     if (ancCode < 10000) {
1442                         ancChg  = ancPdg->Charge();
1443                     } else {
1444                         ancChg  = 0.0;
1445                     }
1446                     Float_t           vxAnc   = ancPart->Vx();
1447                     Float_t           vyAnc   = ancPart->Vy();
1448                     Float_t           vrAnc   = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1449                     fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1450                     ancestor = ancPart->GetFirstMother();           //Get the next ancestor
1451                 }
1452
1453                 //Determine the origin point of the primary particle associated with the hit
1454                 TParticle     *primPart = gAlice->Particle(idprim);
1455                 TParticlePDG  *primPdg  = primPart->GetPDG();
1456                 Int_t      primCode = primPart->GetPdgCode();
1457                 Double_t           primChg;
1458                 if (primCode < 10000) {
1459                     primChg  = primPdg->Charge();
1460                 } else {
1461                     primChg  = 0.0;
1462                 }
1463                 Float_t    vxPrim = primPart->Vx();
1464                 Float_t    vyPrim = primPart->Vy();
1465                 Float_t    vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1466                 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1467             } // Hit Loop
1468         } //if (pEMCAL)
1469     } // Track Loop
1470 }
1471
1472 Int_t AliEMCALJetFinder
1473 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1474
1475     Int_t flag    = 0;
1476     Int_t parton  = 0;
1477     Int_t neutral = 0;
1478
1479     if (charge == 0) neutral = 1;
1480
1481     if (TMath::Abs(code) <= 6   ||
1482         code == 21              ||
1483         code == 92) parton = 1;
1484
1485     //It's not a parton, it's charged and it went through the TPC
1486     if (!parton && !neutral && radius <= 84.0) flag = 1;
1487
1488     return flag;
1489 }
1490
1491
1492
1493 void AliEMCALJetFinder::SaveBackgroundEvent()
1494 {
1495 // Saves the eta-phi lego and the tracklist
1496 //
1497     if (fLegoB) {
1498        fLegoB->Reset();
1499        (*fLegoB) = (*fLegoB) + (*fLego);
1500        if(fDebug)
1501        printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
1502        fLegoB->Integral(), fLego->Integral());
1503     }
1504
1505     if (fPtB)        delete[] fPtB;
1506     if (fEtaB)       delete[] fEtaB;
1507     if (fPhiB)       delete[] fPhiB;
1508     if (fPdgB)       delete[] fPdgB;
1509     if (fTrackListB) delete[] fTrackListB;
1510
1511     fPtB          = new Float_t[fNtS];
1512     fEtaB         = new Float_t[fNtS];
1513     fPhiB         = new Float_t[fNtS];
1514     fPdgB         = new Int_t  [fNtS];
1515     fTrackListB   = new Int_t  [fNtS];
1516
1517     fNtB = 0;
1518
1519     for (Int_t i = 0; i < fNt; i++) {
1520         if (!fTrackList[i]) continue;
1521         fPtB [fNtB]       = fPtT [i];
1522         fEtaB[fNtB]       = fEtaT[i];
1523         fPhiB[fNtB]       = fPhiT[i];
1524         fPdgB[fNtB]       = fPdgT[i];
1525         fTrackListB[fNtB] = 1;
1526         fNtB++;
1527     }
1528     fBackground = 1;
1529     printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
1530 }
1531
1532 void AliEMCALJetFinder::InitFromBackground()
1533 {
1534 //
1535 //
1536     if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
1537
1538     if (fLego) {
1539         fLego->Reset();
1540         (*fLego) = (*fLego) + (*fLegoB);
1541         if(fDebug)
1542             printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
1543                    fLego->Integral(), fLegoB->Integral());
1544     } else {
1545         printf(" => fLego undefined \n");
1546     }
1547 }
1548
1549
1550 void AliEMCALJetFinder::FindTracksInJetCone()
1551 {
1552 //
1553 //  Build list of tracks inside jet cone
1554 //
1555 //  Loop over jets
1556     Int_t njet = Njets();
1557     for (Int_t nj = 0; nj < njet; nj++)
1558     {
1559         Float_t etaj = JetEtaW(nj);
1560         Float_t phij = JetPhiW(nj);
1561         Int_t   nT   = 0;           // number of associated tracks
1562
1563 // Loop over particles in current event
1564         for (Int_t part = 0; part < fNt; part++) {
1565             if (!fTrackList[part]) continue;
1566             Float_t phi      = fPhiT[part];
1567             Float_t eta      = fEtaT[part];
1568             Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1569                                            (phij-phi)*(phij-phi));
1570             if (dr < fConeRadius) {
1571                 fTrackList[part] = nj+2;
1572                 nT++;
1573             } // < ConeRadius ?
1574         } // particle loop
1575
1576 // Same for background event if available
1577         Int_t nTB = 0;
1578         if (fBackground) {
1579             for (Int_t part = 0; part < fNtB; part++) {
1580                 Float_t phi      = fPhiB[part];
1581                 Float_t eta      = fEtaB[part];
1582                 Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1583                                                (phij-phi)*(phij-phi));
1584                 fTrackListB[part] = 1;
1585
1586                 if (dr < fConeRadius) {
1587                     fTrackListB[part] = nj+2;
1588                     nTB++;
1589                 } // < ConeRadius ?
1590             } // particle loop
1591         } // Background available ?
1592
1593         Int_t nT0 = nT + nTB;
1594         printf("Total number of tracks %d\n", nT0);
1595
1596         if (nT0 > 1000) nT0 = 1000;
1597
1598         Float_t* ptT  = new Float_t[nT0];
1599         Float_t* etaT = new Float_t[nT0];
1600         Float_t* phiT = new Float_t[nT0];
1601         Int_t*   pdgT = new Int_t[nT0];
1602
1603         Int_t iT = 0;
1604         Int_t j;
1605
1606         for (Int_t part = 0; part < fNt; part++) {
1607             if (fTrackList[part] == nj+2) {
1608                 Int_t index = 0;
1609                 for (j=iT-1; j>=0; j--) {
1610                     if (fPtT[part] > ptT[j]) {
1611                         index = j+1;
1612                         break;
1613                     }
1614                 }
1615                 for (j=iT-1; j>=index; j--) {
1616                     ptT [j+1]  = ptT [j];
1617                     etaT[j+1]  = etaT[j];
1618                     phiT[j+1]  = phiT[j];
1619                     pdgT[j+1]  = pdgT[j];
1620                 }
1621                 ptT [index] = fPtT [part];
1622                 etaT[index] = fEtaT[part];
1623                 phiT[index] = fPhiT[part];
1624                 pdgT[index] = fPdgT[part];
1625                 iT++;
1626             } // particle associated
1627             if (iT > nT0) break;
1628         } // particle loop
1629
1630         if (fBackground) {
1631             for (Int_t part = 0; part < fNtB; part++) {
1632                 if (fTrackListB[part] == nj+2) {
1633                     Int_t index = 0;
1634                     for (j=iT-1; j>=0; j--) {
1635                         if (fPtB[part] > ptT[j]) {
1636                             index = j+1;
1637
1638                             break;
1639                         }
1640                     }
1641                     for (j=iT-1; j>=index; j--) {
1642                         ptT [j+1]  = ptT [j];
1643                         etaT[j+1]  = etaT[j];
1644                         phiT[j+1]  = phiT[j];
1645                         pdgT[j+1]  = pdgT[j];
1646                     }
1647                     ptT [index] = fPtB [part];
1648                     etaT[index] = fEtaB[part];
1649                     phiT[index] = fPhiB[part];
1650                     pdgT[index] = fPdgB[part];
1651                     iT++;
1652                 } // particle associated
1653                 if (iT > nT0) break;
1654             } // particle loop
1655         } // Background available ?
1656
1657         fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
1658         delete[] ptT;
1659         delete[] etaT;
1660         delete[] phiT;
1661         delete[] pdgT;
1662
1663     } // jet loop loop
1664 }
1665
1666 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1667 {
1668 // Propagates phi angle to EMCAL radius
1669 //
1670   static Float_t b = 0.0, rEMCAL = -1.0;
1671 // Get field in kGS
1672   b = gAlice->Field()->SolenoidField();
1673 // Get EMCAL radius in cm
1674   rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1675   Float_t dPhi = 0.;
1676 //
1677 //
1678 // bending radies
1679 // pt [Gev]
1680 // B  [kG]
1681 //
1682     Float_t rB = 3335.6 * pt / b;  // [cm]  (case of |charge|=1)
1683 //
1684 // check if particle is curling below EMCAL
1685     if (2.*rB < rEMCAL) {
1686         curls = kTRUE;
1687         return dPhi;
1688     }
1689 //
1690 // if not calculate delta phi
1691     Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1692     dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1693     dPhi = -TMath::Sign(dPhi, charge);
1694 //
1695     return dPhi;
1696 }
1697
1698 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1699 {
1700 // dummy for hbook calls
1701     ;
1702 }
1703
1704 void AliEMCALJetFinder::DrawLego(Char_t *opt)
1705 {fLego->Draw(opt);}
1706
1707 void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt)
1708 {fhLegoEMCAL->Draw(opt);}
1709
1710 void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1711 {
1712   static TPaveText *varLabel=0;
1713   if(!fC1) {
1714     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1715   }
1716   fC1->Clear();
1717   fC1->Divide(2,2);
1718   fC1->cd(1);
1719   gPad->SetLogy(1);
1720   fhCellEt->Draw();
1721
1722   fC1->cd(2);
1723   gPad->SetLogy(1);
1724   fhCellEMCALEt->Draw();
1725
1726   fC1->cd(3);
1727   gPad->SetLogy(1);
1728   fhTrackPt->Draw();
1729   fhTrackPtBcut->SetLineColor(2);
1730   fhTrackPtBcut->Draw("same");
1731
1732   fC1->cd(4);
1733   if(!varLabel) {
1734     PrintParameters(1);
1735     varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1736     varLabel->SetTextAlign(12);
1737     varLabel->SetFillColor(19); // see TAttFill
1738     TString tmp(GetTitle());
1739     varLabel->ReadFile(GetFileNameForParameters());
1740   }
1741   varLabel->Draw();
1742   fC1->Update();
1743   if(mode) { // for saving picture to the file
1744     TString stmp(GetFileNameForParameters());
1745     stmp.ReplaceAll("_Par.txt",".ps");
1746     fC1->Print(stmp.Data());
1747   }
1748 }
1749
1750 void AliEMCALJetFinder::PrintParameters(Int_t mode)
1751 {
1752   FILE *file=0;
1753   if(mode==0) file = stdout; // output to terminal
1754   else {
1755     file = fopen(GetFileNameForParameters(),"w");
1756     if(file==0) file = stdout;
1757   }
1758   fprintf(file,"====   Filling lego   ==== \n");
1759   fprintf(file,"Smearing          %6i  ", fSmear);
1760   fprintf(file,"Efficiency        %6i\n", fEffic);
1761   fprintf(file,"Hadr.Correct.     %6i  ", fHCorrection);
1762   fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1763   fprintf(file,"====  Jet finding     ==== \n");
1764   fprintf(file,"Cone radius       %6.2f  ", fConeRadius);
1765   fprintf(file,"Seed E_{T}           %6.1f\n", fEtSeed);
1766   fprintf(file,"Min E_{T} of cell    %6.1f  ", fMinCellEt);
1767   fprintf(file,"Min E_{T} of jet     %6.1f\n", fMinJetEt);
1768   if(fMode) {
1769     fprintf(file,"====  Bg subtraction     ==== \n");
1770     fprintf(file,"BG subtraction  %6i  ", fMode);
1771     fprintf(file,"Min cone move   %6.2f\n", fMinMove);
1772     fprintf(file,"Max cone move   %6.2f  ", fMaxMove);
1773     fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1774   } else
1775   fprintf(file,"==== No Bg subtraction     ==== \n");
1776   if(file != stdout) fclose(file);
1777 }
1778
1779 void AliEMCALJetFinder::DrawLegos()
1780 {
1781   if(!fC1) {
1782     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1783   }
1784   fC1->Clear();
1785   fC1->Divide(2,2);
1786   gStyle->SetOptStat(111111);
1787
1788   Int_t nent1, nent2, nent3, nent4;
1789   Double_t int1, int2, int3, int4;
1790   nent1 = (Int_t)fLego->GetEntries();
1791   int1  = fLego->Integral();
1792   fC1->cd(1);
1793   if(int1) fLego->Draw("lego");
1794
1795   nent2 = (Int_t)fhLegoTracks->GetEntries();
1796   int2  = fhLegoTracks->Integral();
1797   fC1->cd(2);
1798   if(int2) fhLegoTracks->Draw("lego");
1799
1800   nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1801   int3  = fhLegoEMCAL->Integral();
1802   fC1->cd(3);
1803   if(int3) fhLegoEMCAL->Draw("lego");
1804
1805   nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1806   int4  = fhLegoHadrCorr->Integral();
1807   fC1->cd(4);
1808   if(int4) fhLegoHadrCorr->Draw("lego");
1809
1810   // just for checking
1811   printf(" Integrals \n");
1812   printf("lego   %10.3f \ntrack  %10.3f \nhits   %10.3f \nHCorr   %10.3f\n--      %10.3f(must be 0)\n",
1813   int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1814 }
1815
1816 const Char_t* AliEMCALJetFinder::GetFileNameForParameters(Char_t* dir)
1817 {
1818   static TString tmp;
1819   if(strlen(dir)) tmp = dir;
1820   tmp += GetTitle();
1821   tmp += "_Par.txt";
1822   return tmp.Data();
1823 }
1824
1825 void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1826 { // See FillFromTracks() - npart must be positive
1827     if (fTrackList) delete[] fTrackList;
1828     if (fPtT)       delete[] fPtT;
1829     if (fEtaT)      delete[] fEtaT;
1830     if (fPhiT)      delete[] fPhiT;
1831     if (fPdgT)      delete[] fPdgT;
1832
1833     if(npart>0) {
1834        fTrackList = new Int_t  [npart];
1835        fPtT       = new Float_t[npart];
1836        fEtaT      = new Float_t[npart];
1837        fPhiT      = new Float_t[npart];
1838        fPdgT      = new Int_t[npart];
1839     } else {
1840        printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1841     }
1842 }
1843
1844 Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1845 {
1846   Int_t absPdg = TMath::Abs(pdg);
1847   if(absPdg<=6) return kTRUE; // quarks
1848   if(pdg == 21) return kTRUE; // gluon
1849   if(pdg == 92) return kTRUE; // string
1850
1851   // see p.51 of Pythia Manual
1852   // Not include diquarks with c and b quark - 4-mar-2002
1853   //                 ud_0,sd_0,su_0; dd_1,ud_1,uu_1;  sd_1,su_1,ss_1
1854   static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203,  3103,3203,3303};
1855   for(Int_t i=0; i<9; i++) if(absPdg == diquark[i])  return kTRUE; // diquarks
1856
1857   return kFALSE;
1858 }
1859
1860 void AliEMCALJetFinder::FindChargedJet()
1861 {
1862 //
1863 // Find jet from charged particle information only
1864 //
1865
1866 //
1867 //  Look for seeds
1868 //
1869     Int_t njets = 0;
1870     Int_t part  = 0;
1871     Int_t nseed = 0;
1872
1873 //
1874 //
1875     ResetJets();
1876
1877 //
1878     for (part = 0; part < fNt; part++) {
1879         if (!fTrackList[part]) continue;
1880         if (fPtT[part] > fEtSeed) nseed++;
1881     }
1882     printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
1883     Int_t* iSeeds = new Int_t[nseed];
1884     nseed = 0;
1885
1886     for (part = 0; part < fNt; part++) {
1887         if (!fTrackList[part]) continue;
1888         if (fPtT[part] > fEtSeed) iSeeds[nseed++] =  part;
1889     }
1890
1891 //
1892 // Loop over seeds
1893 //
1894     Int_t seed = 0;
1895     Float_t pt;
1896
1897     while(1){
1898 //
1899 // Find seed with highest pt
1900 //
1901         Float_t ptmax = -1.;
1902         Int_t   index = -1;
1903         Int_t   jndex = -1;
1904         for (seed = 0; seed < nseed; seed++) {
1905             if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
1906                 ptmax = pt;
1907                 index = seed;
1908             } // ptmax ?
1909         } // seeds
1910         if (ptmax < 0.) break;
1911         jndex = iSeeds[index];
1912 //
1913 // Remove from the list
1914         iSeeds[index] = -1;
1915         printf("\n Next Seed %d %f", jndex, ptmax);
1916 //
1917 // Find tracks in cone around seed
1918 //
1919         Float_t phiSeed = fPhiT[jndex];
1920         Float_t etaSeed = fEtaT[jndex];
1921         Float_t eT = 0.;
1922         Float_t pxJ = 0.;
1923         Float_t pyJ = 0.;
1924         Float_t pzJ = 0.;
1925
1926         for (part = 0; part < fNt; part++) {
1927             if (!fTrackList[part]) continue;
1928             Float_t deta = fEtaT[part] - etaSeed;
1929             Float_t dphi = fPhiT[part] - phiSeed;
1930             Float_t dR   = TMath::Sqrt(deta * deta + dphi * dphi);
1931             if (dR < fConeRadius) {
1932                 eT += fPtT[part];
1933                 Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
1934                 Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
1935                 Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
1936                 Float_t pz = fPtT[part] / TMath::Tan(theta);
1937                 pxJ += px;
1938                 pyJ += py;
1939                 pzJ += pz;
1940                 //
1941                 // if seed, remove it
1942                 //
1943                 for (seed = 0; seed < nseed; seed++) {
1944                     if (part == iSeeds[seed]) iSeeds[seed] = -1;
1945                 } // seed ?
1946             } // < cone radius
1947         } // particle loop
1948
1949 //
1950 //      Estimate of jet direction
1951         Float_t phiJ   = TMath::ATan2(pyJ, pxJ);
1952         Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
1953         Float_t etaJ   = TMath::Log(TMath::Tan(thetaJ / 2.));
1954 //      Float_t ptJ    = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
1955
1956 //
1957 //      Sum up all energy
1958 //
1959         Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
1960         Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
1961         Int_t dIphi = Int_t(fConeRadius / fDphi);
1962         Int_t dIeta = Int_t(fConeRadius / fDeta);
1963         Int_t iPhi, iEta;
1964         Float_t sumE = 0;
1965         for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
1966             for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
1967                 if (iPhi < 0 || iEta < 0) continue;
1968                 Float_t dPhi = fPhiMin + iPhi * fDphi;
1969                 Float_t dEta = fEtaMin + iEta * fDeta;
1970                 if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
1971                 sumE += fLego->GetBinContent(iEta, iPhi);
1972             } // eta
1973         } // phi
1974 //
1975 //
1976 //
1977         fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
1978         FindTracksInJetCone();
1979         printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
1980         printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
1981         printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
1982     } // while(1)
1983     EMCALJETS.njet = njets;
1984     if (fWrite) WriteJets();
1985     fEvent++;
1986 }