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