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