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