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