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