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