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