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