EMCAL/AliEMCALJetFinder.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /*
  18 $Log$
  19 Revision 1.4  2002/01/21 11:03:21  morsch
  20 Phi propagation introduced in FillFromTracks.
  21
  22 Revision 1.3  2002/01/18 05:07:56  morsch
  23 - hadronic correction
  24 - filling of digits
  25 - track selection upon EMCAL information
  26
  27 */
  28
  29 //*-- Author: Andreas Morsch (CERN)
  30 //*           J.L. Klay      (LBL)
  31 // From root ...
  32 #include <TClonesArray.h>
  33 #include <TTree.h>
  34 #include <TBranchElement.h>
  35 #include <TFile.h>
  36 #include <TH2.h>
  37 #include <TAxis.h>
  38 #include <TParticle.h>
  39 #include <TParticlePDG.h>
  40
  41 // From AliRoot ...
  42 #include "AliEMCALJetFinder.h"
  43 #include "AliEMCALFast.h"
  44 #include "AliEMCALGeometry.h"
  45 #include "AliEMCALHit.h"
  46 #include "AliEMCALDigit.h"
  47 #include "AliEMCALDigitizer.h"
  48 #include "AliEMCALHadronCorrection.h"
  49 #include "Ecommon.h"
  50 #include "AliRun.h"
  51 #include "AliMagF.h"
  52 #include "AliMagFCM.h"
  53 #include "AliEMCAL.h"
  54 #include "AliHeader.h"
  55 #include "AliPDG.h"
  56
  57 ClassImp(AliEMCALJetFinder)
  58
  59 //____________________________________________________________________________
  60 AliEMCALJetFinder::AliEMCALJetFinder()
  61 {
  62 // Default constructor
  63     fJets             = 0;
  64     fNjets            = 0;
  65     fLego             = 0;
  66     fTrackList        = 0;
  67     fPtT              = 0;
  68     fEtaT             = 0;
  69     fPhiT             = 0;
  70     fHadronCorrector  = 0;
  71 }
  72
  73 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
  74     : TTask(name, title)
  75 {
  76 // Constructor
  77     fJets  = new TClonesArray("AliEMCALJet",10000);
  78     fNjets = 0;
  79     for (Int_t i = 0; i < 30000; i++)
  80     {
  81         fEtCell[i]  = 0.;
  82         fEtaCell[i] = 0.;
  83         fPhiCell[i] = 0.;
  84     }
  85     fLego = 0;
  86     fTrackList = 0;
  87     fTrackList = 0;
  88     fPtT       = 0;
  89     fEtaT      = 0;
  90     fPhiT      = 0;
  91     fHadronCorrector = 0;
  92     SetPtCut();
  93     SetMomentumSmearing();
  94     SetEfficiencySim();
  95     SetDebug();
  96     SetHadronCorrection();
  97     SetSamplingFraction();
  98     SetIncludeK0andN();
  99 }
 100
 101
 102
 103
 104 //____________________________________________________________________________
 105 AliEMCALJetFinder::~AliEMCALJetFinder()
 106 {
 107 // Destructor
 108 //
 109     if (fJets){
 110         fJets->Delete();
 111         delete fJets;
 112     }
 113 }
 114
 115
 116
 117
 118 #ifndef WIN32
 119 # define jet_finder_ua1 jet_finder_ua1_
 120 # define hf1 hf1_
 121 # define type_of_call
 122
 123 #else
 124 # define jet_finder_ua1 J
 125 # define hf1 HF1
 126 # define type_of_call _stdcall
 127 #endif
 128
 129 extern "C" void type_of_call
 130 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
 131                Float_t  etc[30000],  Float_t etac[30000],
 132                Float_t  phic[30000],
 133                Float_t& min_move, Float_t& max_move, Int_t& mode,
 134                Float_t& prec_bg,  Int_t& ierror);
 135
 136 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
 137
 138
 139
 140
 141
 142 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
 143                              Float_t etac[30000], Float_t phic[30000],
 144                              Float_t min_move, Float_t max_move, Int_t mode,
 145                              Float_t prec_bg,  Int_t   ierror)
 146 {
 147 // Wrapper for fortran coded jet finder
 148 // Full argument list
 149     jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
 150                    min_move, max_move, mode, prec_bg, ierror);
 151     // Write result to output
 152     WriteJets();
 153 }
 154
 155 void AliEMCALJetFinder::Find()
 156 {
 157 // Wrapper for fortran coded jet finder using member data for
 158 // argument list
 159
 160     Float_t min_move = 0;
 161     Float_t max_move = 0;
 162     Int_t   mode     = 0;
 163     Float_t prec_bg  = 0.;
 164     Int_t   ierror   = 0;
 165
 166     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
 167                    min_move, max_move, mode, prec_bg, ierror);
 168     // Write result to output
 169     Int_t njet = Njets();
 170
 171     for (Int_t nj=0; nj<njet; nj++)
 172     {
 173
 174
 175         fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
 176                                     JetPhiW(nj),
 177                                     JetEtaW(nj));
 178     }
 179     FindTracksInJetCone();
 180     WriteJets();
 181 }
 182
 183
 184 Int_t AliEMCALJetFinder::Njets()
 185 {
 186 // Get number of reconstructed jets
 187     return EMCALJETS.njet;
 188 }
 189
 190 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
 191 {
 192 // Get reconstructed jet energy
 193     return EMCALJETS.etj[i];
 194 }
 195
 196 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
 197 {
 198 // Get reconstructed jet phi from leading particle
 199     return EMCALJETS.phij[0][i];
 200 }
 201
 202 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
 203 {
 204 // Get reconstructed jet phi from weighting
 205     return EMCALJETS.phij[1][i];
 206 }
 207
 208 Float_t  AliEMCALJetFinder::JetEtaL(Int_t i)
 209 {
 210 // Get reconstructed jet eta from leading particles
 211     return EMCALJETS.etaj[0][i];
 212 }
 213
 214
 215 Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)
 216 {
 217 // Get reconstructed jet eta from weighting
 218     return EMCALJETS.etaj[1][i];
 219 }
 220
 221 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
 222 {
 223 // Set grid cell size
 224     EMCALCELLGEO.etaCellSize = eta;
 225     EMCALCELLGEO.phiCellSize = phi;
 226     fDeta = eta;
 227     fDphi = phi;
 228 }
 229
 230 void AliEMCALJetFinder::SetConeRadius(Float_t par)
 231 {
 232 // Set jet cone radius
 233     EMCALJETPARAM.coneRad = par;
 234     fConeRadius = par;
 235 }
 236
 237 void AliEMCALJetFinder::SetEtSeed(Float_t par)
 238 {
 239 // Set et cut for seeds
 240     EMCALJETPARAM.etSeed = par;
 241     fEtSeed = par;
 242 }
 243
 244 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
 245 {
 246 // Set minimum jet et
 247     EMCALJETPARAM.ejMin = par;
 248     fMinJetEt = par;
 249 }
 250
 251 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
 252 {
 253 // Set et cut per cell
 254     EMCALJETPARAM.etMin = par;
 255     fMinCellEt = par;
 256 }
 257
 258 void AliEMCALJetFinder::SetPtCut(Float_t par)
 259 {
 260 // Set pT cut on charged tracks
 261     fPtCut = par;
 262 }
 263
 264
 265 void AliEMCALJetFinder::Test()
 266 {
 267 //
 268 // Test the finder call
 269 //
 270     const Int_t nmax = 30000;
 271     Int_t ncell      = 10;
 272     Int_t ncell_tot  = 100;
 273
 274     Float_t etc[nmax];
 275     Float_t etac[nmax];
 276     Float_t phic[nmax];
 277     Float_t min_move = 0;
 278     Float_t max_move = 0;
 279     Int_t   mode = 0;
 280     Float_t prec_bg = 0;
 281     Int_t   ierror = 0;
 282
 283
 284     Find(ncell, ncell_tot, etc, etac, phic,
 285          min_move, max_move, mode, prec_bg, ierror);
 286
 287 }
 288
 289 //
 290 //  I/O
 291 //
 292
 293 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
 294 {
 295     //
 296     // Add a jet
 297     //
 298     TClonesArray &lrawcl = *fJets;
 299     new(lrawcl[fNjets++]) AliEMCALJet(jet);
 300 }
 301
 302 void AliEMCALJetFinder::ResetJets()
 303 {
 304     //
 305     // Reset Jet List
 306     //
 307     fJets->Clear();
 308     fNjets = 0;
 309 }
 310
 311 void AliEMCALJetFinder::WriteJets()
 312 {
 313 //
 314 // Add all jets to the list
 315 //
 316     const Int_t kBufferSize = 4000;
 317     TTree *pK = gAlice->TreeK();
 318     const char* file = (pK->GetCurrentFile())->GetName();
 319 // I/O
 320     AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
 321
 322     if (fDebug > 1)
 323     printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
 324
 325     if (fJets && gAlice->TreeR()) {
 326         pEMCAL->MakeBranchInTree(gAlice->TreeR(),
 327                                  "EMCALJets",
 328                                  &fJets,
 329                                  kBufferSize,
 330                                  file);
 331     }
 332     Int_t njet = Njets();
 333     for (Int_t nj = 0; nj < njet; nj++)
 334     {
 335         AddJet(*fJetT[nj]);
 336         delete fJetT[nj];
 337     }
 338
 339     Int_t nev = gAlice->GetHeader()->GetEvent();
 340     gAlice->TreeR()->Fill();
 341     char hname[30];
 342     sprintf(hname,"TreeR%d", nev);
 343     gAlice->TreeR()->Write(hname);
 344     gAlice->TreeR()->Reset();
 345     ResetJets();
 346 }
 347
 348 void AliEMCALJetFinder::BookLego()
 349 {
 350 //
 351 //  Book histo for discretisation
 352 //
 353 //
 354 //  Get geometry parameters from
 355     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 356     AliEMCALGeometry* geom =
 357         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 358     fNbinEta = geom->GetNZ();
 359     fNbinPhi = geom->GetNPhi();
 360     const Float_t  phiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
 361     const Float_t  phiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
 362     fDphi   = (phiMax-phiMin)/fNbinEta;
 363     fDeta   = 1.4/fNbinEta;
 364     fNtot   = fNbinPhi*fNbinEta;
 365
 366 //
 367     fLego = new TH2F("legoH","eta-phi",
 368                            fNbinEta, -0.7,  0.7,
 369                            fNbinPhi, phiMin, phiMax);
 370 }
 371
 372 void AliEMCALJetFinder::DumpLego()
 373 {
 374 //
 375 // Dump lego histo into array
 376 //
 377     fNcell = 0;
 378     for (Int_t i = 1; i < fNbinEta; i++) {
 379         for (Int_t j = 1; j < fNbinPhi; j++) {
 380             Float_t e    = fLego->GetBinContent(i,j);
 381             TAxis* Xaxis = fLego->GetXaxis();
 382             TAxis* Yaxis = fLego->GetYaxis();
 383             Float_t eta  = Xaxis->GetBinCenter(i);
 384             Float_t phi  = Yaxis->GetBinCenter(j);
 385             fEtCell[fNcell]  = e;
 386             fEtaCell[fNcell] = eta;
 387             fPhiCell[fNcell] = phi;
 388             fNcell++;
 389         } // phi
 390     } // eta
 391     fNcell--;
 392 }
 393
 394 void AliEMCALJetFinder::ResetMap()
 395 {
 396 //
 397 // Reset eta-phi array
 398
 399     for (Int_t i=0; i<30000; i++)
 400     {
 401         fEtCell[i]  = 0.;
 402         fEtaCell[i] = 0.;
 403         fPhiCell[i] = 0.;
 404     }
 405 }
 406
 407
 408 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 409 {
 410 //
 411 // Fill Cells with track information
 412 //
 413 //
 414     ResetMap();
 415
 416     if (!fLego) BookLego();
 417 // Reset
 418     if (flag == 0) fLego->Reset();
 419 //
 420 // Access particle information
 421     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 422 // Create track list
 423 //
 424 // 0: not selected
 425 // 1: selected for jet finding
 426 // 2: inside jet
 427 // ....
 428     if (fTrackList) delete[] fTrackList;
 429     if (fPtT)       delete[] fPtT;
 430     if (fEtaT)      delete[] fEtaT;
 431     if (fPhiT)      delete[] fPhiT;
 432
 433     fTrackList = new Int_t  [npart];
 434     fPtT       = new Float_t[npart];
 435     fEtaT      = new Float_t[npart];
 436     fPhiT      = new Float_t[npart];
 437     fNt        = npart;
 438     for (Int_t part = 2; part < npart; part++) {
 439         TParticle *MPart = gAlice->Particle(part);
 440         Int_t mpart   = MPart->GetPdgCode();
 441         Int_t child1  = MPart->GetFirstDaughter();
 442         Float_t pT    = MPart->Pt();
 443         Float_t p     = MPart->P();
 444         Float_t phi   = MPart->Phi();
 445         Float_t eta   = MPart->Eta();
 446
 447         if (fDebug > 0) {
 448             if (part == 6 || part == 7)
 449             {
 450                 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f",
 451                        part-5, pT, eta, phi);
 452             }
 453
 454 //          if (mpart == 21)
 455
 456 //              printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
 457 //                     part, mpart, pT, eta, phi);
 458         }
 459
 460
 461         fTrackList[part] = 0;
 462         fPtT[part]       = pT;
 463         fEtaT[part]      = eta;
 464         fPhiT[part]      = phi;
 465
 466         if (part < 2) continue;
 467         if (pT == 0 || pT < fPtCut) continue;
 468         TParticlePDG* pdgP = 0;
 469 // charged or neutral
 470         if (ich == 0) {
 471             pdgP = MPart->GetPDG();
 472             if (pdgP->Charge() == 0) {
 473                 if (!fK0N) {
 474                     continue;
 475                 } else {
 476                     if (mpart != kNeutron    &&
 477                         mpart != kNeutronBar &&
 478                         mpart != kK0Long) continue;
 479                 }
 480             }
 481         }
 482 // skip partons
 483         if (TMath::Abs(mpart) <= 6         ||
 484             mpart == 21                    ||
 485             mpart == 92) continue;
 486 // acceptance cut
 487         if (TMath::Abs(eta) > 0.7)         continue;
 488         if (phi*180./TMath::Pi() > 120.)   continue;
 489 // final state only
 490         if (child1 >= 0 && child1 < npart) continue;
 491 //
 492         if (fDebug > 1)
 493         printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
 494         part, mpart, child1, eta, phi, pT);
 495 //
 496 //
 497 // phi propagation
 498
 499         Bool_t curls = kFALSE;
 500         Float_t dphi = PropagatePhi(pT, pdgP->Charge(), curls);
 501         if (curls) continue;
 502         phi += dphi;
 503 //
 504 // Momentum smearing goes here ...
 505 //
 506         if (fSmear) {
 507             p = AliEMCALFast::SmearMomentum(1,p);
 508         }
 509 //
 510 // Tracking Efficiency and TPC acceptance goes here ...
 511         if (fEffic) {
 512             Float_t eff =  AliEMCALFast::Efficiency(1,p);
 513             if (AliEMCALFast::RandomReject(eff)) continue;
 514         }
 515 //
 516 // Correction of Hadronic Energy goes here ...
 517 //
 518         if (fHCorrection) {
 519             if (!fHadronCorrector)
 520                 Fatal("AliEMCALJetFinder",
 521                     "Hadronic energy correction required but not defined !");
 522             Float_t dpH = fHadronCorrector->GetEnergy(p, eta, 7);
 523             p -= dpH;
 524         }
 525 //
 526 //  More to do ? Just think about it !
 527 //
 528         fTrackList[part] = 1;
 529 //
 530         fLego->Fill(eta, phi, p);
 531     } // primary loop
 532     DumpLego();
 533 }
 534
 535 void AliEMCALJetFinder::FillFromHits(Int_t flag)
 536 {
 537 //
 538 // Fill Cells with hit information
 539 //
 540 //
 541     ResetMap();
 542
 543     if (!fLego) BookLego();
 544     if (flag == 0) fLego->Reset();
 545
 546 //
 547 // Access hit information
 548     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 549
 550     TTree *treeH = gAlice->TreeH();
 551     Int_t ntracks = (Int_t) treeH->GetEntries();
 552 //
 553 //   Loop over tracks
 554 //
 555     Int_t nbytes = 0;
 556
 557
 558     for (Int_t track=0; track<ntracks;track++) {
 559         gAlice->ResetHits();
 560         nbytes += treeH->GetEvent(track);
 561 //
 562 //  Loop over hits
 563 //
 564         for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
 565             mHit;
 566             mHit=(AliEMCALHit*) pEMCAL->NextHit())
 567         {
 568             Float_t x      =    mHit->X();         // x-pos of hit
 569             Float_t y      =    mHit->Y();         // y-pos
 570             Float_t z      =    mHit->Z();         // z-pos
 571             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
 572
 573             Float_t r      =    TMath::Sqrt(x*x+y*y);
 574             Float_t theta  =    TMath::ATan2(r,z);
 575             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
 576             Float_t phi    =    TMath::ATan2(y,x);
 577             fLego->Fill(eta, phi, fSamplingF*eloss);
 578         } // Hit Loop
 579     } // Track Loop
 580     DumpLego();
 581 }
 582
 583 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
 584 {
 585 //
 586 // Fill Cells with digit information
 587 //
 588 //
 589     ResetMap();
 590
 591     if (!fLego) BookLego();
 592     if (flag == 0) fLego->Reset();
 593     Int_t nbytes;
 594
 595
 596 //
 597 //  Connect digits
 598 //
 599     TClonesArray* digs      = new TClonesArray("AliEMCALDigit", 10000);
 600     TTree *treeD = gAlice->TreeD();
 601     TBranchElement* branchDg = (TBranchElement*)
 602         treeD->GetBranch("EMCAL");
 603
 604     if (!branchDg) Fatal("AliEMCALJetFinder",
 605                          "Reading digits requested but no digits in file !");
 606
 607     branchDg->SetAddress(&digs);
 608     Int_t nent = (Int_t) branchDg->GetEntries();
 609 //
 610 //  Connect digitizer
 611 //
 612     AliEMCALDigitizer* digr = new AliEMCALDigitizer();
 613     TBranchElement* branchDr = (TBranchElement*)
 614         treeD->GetBranch("AliEMCALDigitizer");
 615     branchDr->SetAddress(&digr);
 616 //
 617 //
 618     nbytes += branchDg->GetEntry(0);
 619     nbytes += branchDr->GetEntry(0);
 620 //
 621 //  Get digitizer parameters
 622     Float_t towerADCped  = digr->GetTowerpedestal();
 623     Float_t towerADCcha  = digr->GetTowerchannel();
 624     Float_t preshoADCped = digr->GetPreShopedestal();
 625     Float_t preshoADCcha = digr->GetPreShochannel();
 626
 627     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 628     AliEMCALGeometry* geom =
 629         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 630
 631     if (fDebug) {
 632         Int_t   ndig = digs->GetEntries();
 633         printf("\n Number of Digits: %d %d\n", ndig, nent);
 634         printf("\n Parameters: %f %f %f %f\n",
 635                towerADCped, towerADCcha, preshoADCped, preshoADCcha );
 636         printf("\n Geometry: %d %d\n", geom->GetNEta(), geom->GetNPhi());
 637     }
 638
 639 //
 640 //  Loop over digits
 641     AliEMCALDigit* sdg;
 642     TIter next(digs);
 643     while ((sdg = (AliEMCALDigit*)(next())))
 644     {
 645         Double_t pedestal;
 646         Double_t channel;
 647         if (sdg->GetId() > (geom->GetNZ() * geom->GetNPhi()))
 648         {
 649             pedestal = preshoADCped;
 650             channel  = preshoADCcha;
 651         } else {
 652             pedestal = towerADCped;
 653             channel  = towerADCcha;
 654         }
 655
 656         Float_t eta = sdg->GetEta();
 657         Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
 658         Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
 659
 660         if (fDebug > 1) printf("\n Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
 661                            eta, phi, amp, sdg->GetAmp());
 662
 663         fLego->Fill(eta, phi, fSamplingF*amp);
 664     } // digit loop
 665 //
 666 //  Dump lego hist
 667     DumpLego();
 668 }
 669
 670
 671 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
 672 {
 673 //
 674 // Fill Cells with hit information
 675 //
 676 //
 677     ResetMap();
 678
 679     if (!fLego) BookLego();
 680 // Reset
 681     if (flag == 0) fLego->Reset();
 682
 683 //Flag charged tracks passing through TPC which
 684 //are associated to EMCAL Hits
 685     BuildTrackFlagTable();
 686
 687 //
 688 // Access particle information
 689     TTree *treeH = gAlice->TreeH();
 690     Int_t ntracks = (Int_t) treeH->GetEntries();
 691
 692     if (fPtT)       delete[] fPtT;
 693     if (fEtaT)      delete[] fEtaT;
 694     if (fPhiT)      delete[] fPhiT;
 695
 696     fPtT       = new Float_t[ntracks];
 697     fEtaT      = new Float_t[ntracks];
 698     fPhiT      = new Float_t[ntracks];
 699     fNt        = ntracks;
 700
 701     for (Int_t track = 0; track < ntracks; track++) {
 702           TParticle *MPart = gAlice->Particle(track);
 703           Float_t pT    = MPart->Pt();
 704           Float_t phi   = MPart->Phi();
 705           Float_t eta   = MPart->Eta();
 706
 707         if(fTrackList[track]) {
 708           fPtT[track]       = pT;
 709           fEtaT[track]      = eta;
 710           fPhiT[track]      = phi;
 711
 712           if (track < 2) continue;      //Colliding particles?
 713           if (pT == 0 || pT < fPtCut) continue;
 714           fLego->Fill(eta, phi, pT);
 715         }
 716       } // track loop
 717     DumpLego();
 718 }
 719
 720 void AliEMCALJetFinder::BuildTrackFlagTable() {
 721
 722 // Method to generate a lookup table for TreeK particles
 723 // which are linked to hits in the EMCAL
 724 //
 725 // --Author: J.L. Klay
 726 //
 727     printf("\n Building track flag table.\n");
 728
 729
 730 // Access hit information
 731     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 732
 733     TTree *TK = gAlice->TreeK();                // Get the number of entries in the kine tree
 734     Int_t nKTrks = (Int_t) TK->GetEntries();    // (Number of particles created somewhere)
 735
 736     if(fTrackList) delete[] fTrackList;         //Make sure we get rid of the old one
 737     fTrackList = new Int_t[nKTrks];             //before generating a new one
 738
 739     for (Int_t i = 0; i < nKTrks; i++) {        //Initialize members to 0
 740         fTrackList[i] = 0;
 741     }
 742
 743     TTree *treeH = gAlice->TreeH();
 744     Int_t ntracks = (Int_t) treeH->GetEntries();
 745 //
 746 //   Loop over tracks
 747 //
 748     Int_t nbytes = 0;
 749
 750     for (Int_t track=0; track<ntracks;track++) {
 751         gAlice->ResetHits();
 752         nbytes += treeH->GetEvent(track);
 753         if (pEMCAL)  {
 754
 755 //
 756 //  Loop over hits
 757 //
 758             for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
 759                 mHit;
 760                 mHit=(AliEMCALHit*) pEMCAL->NextHit())
 761             {
 762                 Int_t   iTrk   =    mHit->Track();       // track number
 763                 Int_t   idprim =    mHit->GetPrimary();  // primary particle
 764
 765                 //Determine the origin point of this particle - it made a hit in the EMCAL
 766                 TParticle         *trkPart = gAlice->Particle(iTrk);
 767                 TParticlePDG  *trkPdg  = trkPart->GetPDG();
 768                 Int_t      trkCode = trkPart->GetPdgCode();
 769                 Double_t           trkChg;
 770                 if (trkCode < 10000) {          //Big Ions cause problems for
 771                     trkChg  = trkPdg->Charge(); //this function.  Since they aren't
 772                 } else {                                //likely to make it very far, set
 773                     trkChg  = 0.0;                      //their charge to 0 for the Flag test
 774                 }
 775                 Float_t            vxTrk   = trkPart->Vx();
 776                 Float_t            vyTrk   = trkPart->Vy();
 777                 Float_t            vrTrk   = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
 778                 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
 779
 780                 //Loop through the ancestry of the EMCAL entrance particles
 781                 Int_t ancestor = trkPart->GetFirstMother();  //Get track's Mother
 782                 while (ancestor != -1) {
 783                     TParticle     *ancPart = gAlice->Particle(ancestor);  //get particle info on ancestor
 784                     TParticlePDG  *ancPdg  = ancPart->GetPDG();
 785                     Int_t             ancCode = ancPart->GetPdgCode();
 786                     Double_t       ancChg;
 787                     if (ancCode < 10000) {
 788                         ancChg  = ancPdg->Charge();
 789                     } else {
 790                         ancChg  = 0.0;
 791                     }
 792                     Float_t           vxAnc   = ancPart->Vx();
 793                     Float_t           vyAnc   = ancPart->Vy();
 794                     Float_t           vrAnc   = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
 795                     fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
 796                     ancestor = ancPart->GetFirstMother();           //Get the next ancestor
 797                 }
 798
 799                 //Determine the origin point of the primary particle associated with the hit
 800                 TParticle     *primPart = gAlice->Particle(idprim);
 801                 TParticlePDG  *primPdg  = primPart->GetPDG();
 802                 Int_t      primCode = primPart->GetPdgCode();
 803                 Double_t           primChg;
 804                 if (primCode < 10000) {
 805                     primChg  = primPdg->Charge();
 806                 } else {
 807                     primChg  = 0.0;
 808                 }
 809                 Float_t    vxPrim = primPart->Vx();
 810                 Float_t    vyPrim = primPart->Vy();
 811                 Float_t    vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
 812                 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
 813             } // Hit Loop
 814         } //if (pEMCAL)
 815     } // Track Loop
 816 }
 817
 818 Int_t AliEMCALJetFinder
 819 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
 820
 821     Int_t flag = 0;
 822     Int_t parton = 0;
 823     Int_t neutral = 0;
 824
 825     if (charge == 0) neutral = 1;
 826
 827     if (TMath::Abs(code) <= 6   ||
 828         code == 21              ||
 829         code == 92) parton = 1;
 830
 831     //It's not a parton, it's charged and it went through the TPC
 832     if (!parton && !neutral && radius <= 84.0) flag = 1;
 833
 834     return flag;
 835 }
 836
 837 void AliEMCALJetFinder::FindTracksInJetCone()
 838 {
 839 //
 840 //  Build list of tracks inside jet cone
 841 //
 842 //  Loop over jets
 843     Int_t njet = Njets();
 844     for (Int_t nj = 0; nj < njet; nj++)
 845     {
 846         Float_t etaj = JetEtaW(nj);
 847         Float_t phij = JetPhiW(nj);
 848         Int_t   nT   = 0;           // number of associated tracks
 849
 850 // Loop over particles
 851         for (Int_t part = 0; part < fNt; part++) {
 852             if (!fTrackList[part]) continue;
 853             Float_t phi      = fPhiT[part];
 854             Float_t eta      = fEtaT[part];
 855             Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
 856                                            (phij-phi)*(phij-phi));
 857             if (dr < fConeRadius) {
 858                 fTrackList[part] = nj+2;
 859                 nT++;
 860             } // < ConeRadius ?
 861         } // particle loop
 862
 863         if (nT > 50) nT = 50;
 864
 865         Float_t* ptT  = new Float_t[nT];
 866         Float_t* etaT = new Float_t[nT];
 867         Float_t* phiT = new Float_t[nT];
 868         Int_t iT = 0;
 869         Int_t j;
 870
 871         for (Int_t part = 0; part < fNt; part++) {
 872             if (fTrackList[part] == nj+2) {
 873                 Int_t index = 0;
 874                 for (j=iT-1; j>=0; j--) {
 875                     if (fPtT[part] > ptT[j]) {
 876                         index = j+1;
 877                         break;
 878                     }
 879                 }
 880                 for (j=iT-1; j>=index; j--) {
 881                     ptT [j+1]  = ptT [j];
 882                     etaT[j+1]  = etaT[j];
 883                     phiT[j+1]  = phiT[j];
 884                 }
 885                 ptT [index] = fPtT [part];
 886                 etaT[index] = fEtaT[part];
 887                 phiT[index] = fPhiT[part];
 888                 iT++;
 889             } // particle associated
 890         } // particle loop
 891
 892         fJetT[nj]->SetTrackList(nT, ptT, etaT, phiT);
 893         delete[] ptT;
 894         delete[] etaT;
 895         delete[] phiT;
 896     } // jet loop loop
 897 }
 898
 899 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
 900 {
 901 // Propagates phi angle to EMCAL radius
 902 //
 903     Float_t dPhi = 0.;
 904 // Get field
 905     Float_t b =  ((AliMagFCM*) gAlice->Field())->SolenoidField();
 906 // Get EMCAL radius
 907     Float_t rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
 908 //
 909 //
 910 // bending radies
 911     Float_t rB = 3.3356 * pt / b;
 912
 913 //
 914 // check if particle is curling below EMCAL
 915     if (2.*rB < rEMCAL) {
 916         curls = kTRUE;
 917         return dPhi;
 918     }
 919 //
 920 // if not calculate delta phi
 921     Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
 922     dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
 923     dPhi = TMath::Sign(dPhi, charge);
 924 //
 925     return dPhi;
 926
 927 }
 928
 929
 930 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
 931 {
 932 // dummy for hbook calls
 933     ;
 934 }
 935
 936
 937
 938
 939
 940
 941
 942