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