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