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.33  2003/01/10 10:48:19  morsch
  19 SetSamplingFraction() removed from constructor.
  20
  21 Revision 1.32  2003/01/10 10:26:40  morsch
  22 Sampling fraction initialized from geometry class.
  23
  24 Revision 1.31  2003/01/08 17:13:41  schutz
  25 added the HCAL section
  26
  27 Revision 1.30  2002/12/09 16:26:28  morsch
  28 - Nummber of particles per jet increased to 1000
  29 - Warning removed.
  30
  31 Revision 1.29  2002/11/21 17:01:40  alibrary
  32 Removing AliMCProcess and AliMC
  33
  34 Revision 1.28  2002/11/20 14:13:16  morsch
  35 - FindChargedJets() added.
  36 - Destructor corrected.
  37 - Geometry cuts taken from AliEMCALGeometry.
  38
  39 Revision 1.27  2002/11/15 17:39:10  morsch
  40 GetPythiaParticleName removed.
  41
  42 Revision 1.26  2002/10/14 14:55:35  hristov
  43 Merging the VirtualMC branch to the main development branch (HEAD)
  44
  45 Revision 1.20.4.3  2002/10/10 15:07:49  hristov
  46 Updating VirtualMC to v3-09-02
  47
  48 Revision 1.25  2002/09/13 10:24:57  morsch
  49 idem
  50
  51 Revision 1.24  2002/09/13 10:21:13  morsch
  52 No cast to AliMagFCM.
  53
  54 Revision 1.23  2002/06/27 09:24:26  morsch
  55 Uncomment the TH1::AddDirectory statement.
  56
  57 Revision 1.22  2002/05/22 13:48:43  morsch
  58 Pdg code added to track list.
  59
  60 Revision 1.21  2002/04/27 07:43:08  morsch
  61 Calculation of fDphi corrected (Renan Cabrera)
  62
  63 Revision 1.20  2002/03/12 01:06:23  pavlinov
  64 Testin output from generator
  65
  66 Revision 1.19  2002/02/27 00:46:33  pavlinov
  67 Added method FillFromParticles()
  68
  69 Revision 1.18  2002/02/21 08:48:59  morsch
  70 Correction in FillFromHitFlaggedTrack. (Jennifer Klay)
  71
  72 Revision 1.17  2002/02/14 08:52:53  morsch
  73 Major updates by Aleksei Pavlinov:
  74 FillFromPartons, FillFromTracks, jetfinder configuration.
  75
  76 Revision 1.16  2002/02/08 11:43:05  morsch
  77 SetOutputFileName(..) allows to specify an output file to which the
  78 reconstructed jets are written. If not set output goes to input file.
  79 Attention call Init() before processing.
  80
  81 Revision 1.15  2002/02/02 08:37:18  morsch
  82 Formula for rB corrected.
  83
  84 Revision 1.14  2002/02/01 08:55:30  morsch
  85 Fill map with Et and pT.
  86
  87 Revision 1.13  2002/01/31 09:37:36  morsch
  88 Geometry parameters in constructor and call of SetCellSize()
  89
  90 Revision 1.12  2002/01/23 13:40:23  morsch
  91 Fastidious debug print statement removed.
  92
  93 Revision 1.11  2002/01/22 17:25:47  morsch
  94 Some corrections for event mixing and bg event handling.
  95
  96 Revision 1.10  2002/01/22 10:31:44  morsch
  97 Some correction for bg mixing.
  98
  99 Revision 1.9  2002/01/21 16:28:42  morsch
 100 Correct sign of dphi.
 101
 102 Revision 1.8  2002/01/21 16:05:31  morsch
 103 - different phi-bin for hadron correction
 104 - provisions for background mixing.
 105
 106 Revision 1.7  2002/01/21 15:47:18  morsch
 107 Bending radius correctly in cm.
 108
 109 Revision 1.6  2002/01/21 12:53:50  morsch
 110 authors
 111
 112 Revision 1.5  2002/01/21 12:47:47  morsch
 113 Possibility to include K0long and neutrons.
 114
 115 Revision 1.4  2002/01/21 11:03:21  morsch
 116 Phi propagation introduced in FillFromTracks.
 117
 118 Revision 1.3  2002/01/18 05:07:56  morsch
 119 - hadronic correction
 120 - filling of digits
 121 - track selection upon EMCAL information
 122
 123 */
 124
 125 //*-- Authors: Andreas Morsch   (CERN)
 126 //*            J.L. Klay        (LBL)
 127 //*            Aleksei Pavlinov (WSU)
 128
 129 #include <stdio.h>
 130 // From root ...
 131 #include <TArrayF.h>
 132 #include <TAxis.h>
 133 #include <TBranchElement.h>
 134 #include <TCanvas.h>
 135 #include <TClonesArray.h>
 136 #include <TFile.h>
 137 #include <TH1.h>
 138 #include <TH2.h>
 139 #include <TList.h>
 140 #include <TPDGCode.h>
 141 #include <TPad.h>
 142 #include <TParticle.h>
 143 #include <TParticlePDG.h>
 144 #include <TPaveText.h>
 145 #include <TPythia6Calls.h>
 146 #include <TROOT.h>
 147 #include <TStyle.h>
 148 #include <TTree.h>
 149
 150 // From AliRoot ...
 151 #include "AliEMCAL.h"
 152 #include "AliEMCALDigit.h"
 153 #include "AliEMCALDigitizer.h"
 154 #include "AliEMCALFast.h"
 155 #include "AliEMCALGeometry.h"
 156 #include "AliEMCALHadronCorrection.h"
 157 #include "AliEMCALHit.h"
 158 #include "AliEMCALJetFinder.h"
 159 #include "AliEMCALJetMicroDst.h"
 160 #include "AliHeader.h"
 161 #include "AliMagF.h"
 162 #include "AliMagFCM.h"
 163 #include "AliRun.h"
 164
 165 // Interface to FORTRAN
 166 #include "Ecommon.h"
 167
 168
 169 ClassImp(AliEMCALJetFinder)
 170
 171 //____________________________________________________________________________
 172 AliEMCALJetFinder::AliEMCALJetFinder()
 173 {
 174 // Default constructor
 175     fJets             = 0;
 176     fNjets            = 0;
 177     fLego             = 0;
 178     fLegoB            = 0;
 179
 180     fTrackList        = 0;
 181     fPtT              = 0;
 182     fEtaT             = 0;
 183     fPhiT             = 0;
 184     fPdgT             = 0;
 185
 186     fTrackListB       = 0;
 187     fPtB              = 0;
 188     fEtaB             = 0;
 189     fPhiB             = 0;
 190     fPdgB             = 0;
 191
 192     fHCorrection      = 0;
 193     fHadronCorrector  = 0;
 194
 195     fWrite            = 0;
 196     fOutFileName      = 0;
 197     fOutFile          = 0;
 198     fInFile           = 0;
 199     fEvent            = 0;
 200
 201     SetParametersForBgSubtraction();
 202 }
 203
 204 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
 205     : TTask(name, title)
 206 {
 207 // Constructor
 208 // Title is used in method GetFileNameForParameters();
 209 //
 210     fJets  = new TClonesArray("AliEMCALJet",10000);
 211     fNjets = 0;
 212     for (Int_t i = 0; i < 30000; i++)
 213     {
 214         fEtCell[i]  = 0.;
 215         fEtaCell[i] = 0.;
 216         fPhiCell[i] = 0.;
 217     }
 218     fLego       = 0;
 219     fLegoB      = 0;
 220
 221     fTrackList  = 0;
 222     fPtT        = 0;
 223     fEtaT       = 0;
 224     fPhiT       = 0;
 225     fPdgT       = 0;
 226
 227     fTrackListB       = 0;
 228     fPtB        = 0;
 229     fEtaB       = 0;
 230     fPhiB       = 0;
 231     fPdgB       = 0;
 232
 233     fHCorrection      = 0;
 234     fHadronCorrector  = 0;
 235     fBackground       = 0;
 236     fWrite            = 0;
 237     fOutFileName      = 0;
 238     fOutFile          = 0;
 239     fInFile           = 0;
 240     fEvent            = 0;
 241 //
 242     SetPtCut();
 243     SetMomentumSmearing();
 244     SetEfficiencySim();
 245     SetDebug();
 246     SetHadronCorrection();
 247     SetIncludeK0andN();
 248
 249     SetParametersForBgSubtraction();
 250 }
 251
 252 void AliEMCALJetFinder::SetParametersForBgSubtraction
 253 (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
 254 {
 255 // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
 256 // at WSU Linux cluster - 11-feb-2002
 257 // These parameters must be tuned more carefull !!!
 258   SetMode(mode);
 259   SetMinMove(minMove);
 260   SetMaxMove(maxMove);
 261   SetPrecBg(precBg);
 262 }
 263
 264 //____________________________________________________________________________
 265 AliEMCALJetFinder::~AliEMCALJetFinder()
 266 {
 267 // Destructor
 268 //
 269     if (fJets){
 270         fJets->Delete();
 271         delete fJets;
 272     }
 273     delete fLego;
 274     delete fLegoB;
 275     delete fhLegoTracks;
 276     delete fhLegoEMCAL;
 277     delete fhLegoHadrCorr;
 278     delete fhEff;
 279     delete fhCellEt;
 280     delete fhCellEMCALEt;
 281     delete fhTrackPt;
 282     delete fhTrackPtBcut;
 283     delete fhChPartMultInTpc;
 284
 285     delete[] fTrackList;
 286     delete[] fPtT;
 287     delete[] fEtaT;
 288     delete[] fPhiT;
 289     delete[] fPdgT;
 290
 291     delete[] fTrackListB;
 292     delete[] fPtB;
 293     delete[] fEtaB;
 294     delete[] fPhiB;
 295     delete[] fPdgB;
 296 }
 297
 298 #ifndef WIN32
 299 # define jet_finder_ua1 jet_finder_ua1_
 300 # define hf1 hf1_
 301 # define type_of_call
 302
 303 #else
 304 # define jet_finder_ua1 JET_FINDER_UA1
 305 # define hf1 HF1
 306 # define type_of_call _stdcall
 307 #endif
 308
 309 extern "C" void type_of_call
 310 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
 311                Float_t  etc[30000],  Float_t etac[30000],
 312                Float_t  phic[30000],
 313                Float_t& min_move, Float_t& max_move, Int_t& mode,
 314                Float_t& prec_bg,  Int_t& ierror);
 315
 316 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
 317
 318
 319 void AliEMCALJetFinder::Init()
 320 {
 321 //
 322 // Geometry and I/O initialization
 323 //
 324 //
 325 //
 326 //  Get geometry parameters from EMCAL
 327 //
 328 //
 329 //  Geometry
 330     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 331     AliEMCALGeometry* geom =
 332     AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 333
 334     SetSamplingFraction(geom->GetSampling());
 335
 336     fNbinEta = geom->GetNZ();
 337     fNbinPhi = geom->GetNPhi();
 338     fPhiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
 339     fPhiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
 340     fEtaMin  = geom->GetArm1EtaMin();
 341     fEtaMax  = geom->GetArm1EtaMax();
 342     fDphi    = (fPhiMax-fPhiMin)/fNbinPhi;
 343     fDeta    = (fEtaMax-fEtaMin)/fNbinEta;
 344     fNtot    = fNbinPhi*fNbinEta;
 345 //
 346     SetCellSize(fDeta, fDphi);
 347 //
 348 //  I/O
 349     if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
 350 }
 351
 352 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
 353                              Float_t etac[30000], Float_t phic[30000],
 354                              Float_t min_move, Float_t max_move, Int_t mode,
 355                              Float_t prec_bg,  Int_t   ierror)
 356 {
 357 // Wrapper for fortran coded jet finder
 358 // Full argument list
 359     jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
 360                    min_move, max_move, mode, prec_bg, ierror);
 361     // Write result to output
 362     if(fWrite) WriteJets();
 363     fEvent++;
 364 }
 365
 366 void AliEMCALJetFinder::Find()
 367 {
 368 // Wrapper for fortran coded jet finder using member data for
 369 // argument list
 370
 371     Float_t min_move = fMinMove;
 372     Float_t max_move = fMaxMove;
 373     Int_t   mode     = fMode;
 374     Float_t prec_bg  = fPrecBg;
 375     Int_t   ierror;
 376
 377     ResetJets(); // 4-feb-2002 by PAI
 378
 379     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
 380                    min_move, max_move, mode, prec_bg, ierror);
 381     fError = ierror;
 382     // Write result to output
 383     Int_t njet = Njets();
 384
 385     for (Int_t nj=0; nj<njet; nj++)
 386     {
 387
 388         fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
 389                                     JetPhiW(nj),
 390                                     JetEtaW(nj));
 391     }
 392
 393     FindTracksInJetCone();
 394     if(fWrite) WriteJets();
 395     fEvent++;
 396 }
 397
 398 Int_t AliEMCALJetFinder::Njets()
 399 {
 400 // Get number of reconstructed jets
 401     return EMCALJETS.njet;
 402 }
 403
 404 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
 405 {
 406 // Get reconstructed jet energy
 407     return EMCALJETS.etj[i];
 408 }
 409
 410 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
 411 {
 412 // Get reconstructed jet phi from leading particle
 413     return EMCALJETS.phij[0][i];
 414 }
 415
 416 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
 417 {
 418 // Get reconstructed jet phi from weighting
 419     return EMCALJETS.phij[1][i];
 420 }
 421
 422 Float_t  AliEMCALJetFinder::JetEtaL(Int_t i)
 423 {
 424 // Get reconstructed jet eta from leading particles
 425     return EMCALJETS.etaj[0][i];
 426 }
 427
 428
 429 Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)
 430 {
 431 // Get reconstructed jet eta from weighting
 432     return EMCALJETS.etaj[1][i];
 433 }
 434
 435 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
 436 {
 437 // Set grid cell size
 438     EMCALCELLGEO.etaCellSize = eta;
 439     EMCALCELLGEO.phiCellSize = phi;
 440 }
 441
 442 void AliEMCALJetFinder::SetConeRadius(Float_t par)
 443 {
 444 // Set jet cone radius
 445     EMCALJETPARAM.coneRad = par;
 446     fConeRadius = par;
 447 }
 448
 449 void AliEMCALJetFinder::SetEtSeed(Float_t par)
 450 {
 451 // Set et cut for seeds
 452     EMCALJETPARAM.etSeed = par;
 453     fEtSeed = par;
 454 }
 455
 456 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
 457 {
 458 // Set minimum jet et
 459     EMCALJETPARAM.ejMin = par;
 460     fMinJetEt = par;
 461 }
 462
 463 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
 464 {
 465 // Set et cut per cell
 466     EMCALJETPARAM.etMin = par;
 467     fMinCellEt = par;
 468 }
 469
 470 void AliEMCALJetFinder::SetPtCut(Float_t par)
 471 {
 472 // Set pT cut on charged tracks
 473     fPtCut = par;
 474 }
 475
 476
 477 void AliEMCALJetFinder::Test()
 478 {
 479 //
 480 // Test the finder call
 481 //
 482     const Int_t nmax = 30000;
 483     Int_t ncell      = 10;
 484     Int_t ncell_tot  = 100;
 485
 486     Float_t etc[nmax];
 487     Float_t etac[nmax];
 488     Float_t phic[nmax];
 489     Float_t min_move = 0;
 490     Float_t max_move = 0;
 491     Int_t   mode = 0;
 492     Float_t prec_bg = 0;
 493     Int_t   ierror = 0;
 494
 495
 496     Find(ncell, ncell_tot, etc, etac, phic,
 497          min_move, max_move, mode, prec_bg, ierror);
 498
 499 }
 500
 501 //
 502 //  I/O
 503 //
 504
 505 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
 506 {
 507     //
 508     // Add a jet
 509     //
 510     TClonesArray &lrawcl = *fJets;
 511     new(lrawcl[fNjets++]) AliEMCALJet(jet);
 512 }
 513
 514 void AliEMCALJetFinder::ResetJets()
 515 {
 516     //
 517     // Reset Jet List
 518     //
 519     fJets->Clear();
 520     fNjets = 0;
 521 }
 522
 523 void AliEMCALJetFinder::WriteJets()
 524 {
 525 //
 526 // Add all jets to the list
 527 //
 528     const Int_t kBufferSize = 4000;
 529     const char* file = 0;
 530
 531     Int_t njet = Njets();
 532
 533     for (Int_t nj = 0; nj < njet; nj++)
 534     {
 535         AddJet(*fJetT[nj]);
 536         delete fJetT[nj];
 537     }
 538
 539 // I/O
 540     if (!fOutFileName) {
 541 //
 542 // output written to input file
 543 //
 544         AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
 545         TTree* pK = gAlice->TreeK();
 546         file = (pK->GetCurrentFile())->GetName();
 547         if (fDebug > 1)
 548             printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
 549         if (fJets && gAlice->TreeR()) {
 550             pEMCAL->MakeBranchInTree(gAlice->TreeR(),
 551                                      "EMCALJets",
 552                                      &fJets,
 553                                      kBufferSize,
 554                                      file);
 555         }
 556         Int_t nev = gAlice->GetHeader()->GetEvent();
 557         gAlice->TreeR()->Fill();
 558         char hname[30];
 559         sprintf(hname,"TreeR%d", nev);
 560         gAlice->TreeR()->Write(hname);
 561         gAlice->TreeR()->Reset();
 562     } else {
 563 //
 564 // Output written to user specified output file
 565 //
 566         TTree* pK = gAlice->TreeK();
 567         fInFile  = pK->GetCurrentFile();
 568
 569         fOutFile->cd();
 570         char hname[30];
 571         sprintf(hname,"TreeR%d", fEvent);
 572         TTree* treeJ = new TTree(hname, "EMCALJets");
 573         treeJ->Branch("EMCALJets", &fJets, kBufferSize);
 574         treeJ->Fill();
 575         treeJ->Write(hname);
 576         fInFile->cd();
 577     }
 578     ResetJets();
 579 }
 580
 581 void AliEMCALJetFinder::BookLego()
 582 {
 583 //
 584 //  Book histo for discretisation
 585 //
 586
 587 //
 588 //  Don't add histos to the current directory
 589     if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
 590
 591     TH2::AddDirectory(0);
 592     TH1::AddDirectory(0);
 593     gROOT->cd();
 594 //
 595 //  Signal map
 596     fLego = new TH2F("legoH","eta-phi",
 597                            fNbinEta, fEtaMin, fEtaMax,
 598                            fNbinPhi, fPhiMin, fPhiMax);
 599 //
 600 //  Background map
 601     fLegoB = new TH2F("legoB","eta-phi for BG event",
 602                            fNbinEta, fEtaMin, fEtaMax,
 603                            fNbinPhi, fPhiMin, fPhiMax);
 604
 605 //  Tracks map
 606     fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks",
 607     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 608 //  EMCAL map
 609     fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL",
 610     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 611 //  Hadron correction map
 612     fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction",
 613     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 614 //  Hists. for tuning jet finder
 615     fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.);
 616
 617     TArrayF eTmp(1101);
 618     eTmp[0] = 0.0;
 619     for(Int_t i=1; i<=1000; i++)      eTmp[i] = 0.1*i; // step 100 mev
 620     for(Int_t i=1001; i<=1100; i++)   eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV
 621
 622     fhCellEt  = new TH1F("hCellEt","Cell E_{T} from fLego",
 623     eTmp.GetSize()-1, eTmp.GetArray());
 624     fhCellEMCALEt  = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself",
 625     eTmp.GetSize()-1, eTmp.GetArray());
 626     fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ",
 627     eTmp.GetSize()-1, eTmp.GetArray());
 628     fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
 629     eTmp.GetSize()-1, eTmp.GetArray());
 630
 631     fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
 632     "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
 633
 634             //! first canvas for drawing
 635     fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kTRUE);
 636 }
 637
 638 void AliEMCALJetFinder::DumpLego()
 639 {
 640 //
 641 // Dump lego histo into array
 642 //
 643     fNcell = 0;
 644     TAxis* Xaxis = fLego->GetXaxis();
 645     TAxis* Yaxis = fLego->GetYaxis();
 646     //    fhCellEt->Clear();
 647     Float_t e, eH;
 648     for (Int_t i = 1; i <= fNbinEta; i++) {
 649         for (Int_t j = 1; j <= fNbinPhi; j++) {
 650             e = fLego->GetBinContent(i,j);
 651             if (e > 0.0) {
 652               Float_t eta  = Xaxis->GetBinCenter(i);
 653               Float_t phi  = Yaxis->GetBinCenter(j);
 654               fEtCell[fNcell]  = e;
 655               fEtaCell[fNcell] = eta;
 656               fPhiCell[fNcell] = phi;
 657               fNcell++;
 658               fhCellEt->Fill(e);
 659             }
 660             if(fhLegoEMCAL) {
 661               eH = fhLegoEMCAL->GetBinContent(i,j);
 662               if(eH > 0.0) fhCellEMCALEt->Fill(eH);
 663             }
 664         } // phi
 665     } // eta
 666     fNcell--;
 667 }
 668
 669 void AliEMCALJetFinder::ResetMap()
 670 {
 671 //
 672 // Reset eta-phi array
 673
 674     for (Int_t i=0; i<30000; i++)
 675     {
 676         fEtCell[i]  = 0.;
 677         fEtaCell[i] = 0.;
 678         fPhiCell[i] = 0.;
 679     }
 680 }
 681
 682
 683 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 684 {
 685 //
 686 // Fill Cells with track information
 687 //
 688     if (fDebug >= 2)
 689     printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
 690     fNChTpc = 0;
 691
 692     ResetMap();
 693
 694     if (!fLego) BookLego();
 695 // Reset
 696     if (flag == 0) fLego->Reset();
 697 //
 698 // Access particle information
 699     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 700     Int_t ntr   = (gAlice->GetHeader())->GetNtrack();
 701     printf(" : #primary particles %i # tracks %i \n", npart, ntr);
 702
 703 // Create track list
 704 //
 705 // 0: not selected
 706 // 1: selected for jet finding
 707 // 2: inside jet
 708 // ....
 709     if (fTrackList) delete[] fTrackList;
 710     if (fPtT)       delete[] fPtT;
 711     if (fEtaT)      delete[] fEtaT;
 712     if (fPhiT)      delete[] fPhiT;
 713     if (fPdgT)      delete[] fPdgT;
 714
 715     fTrackList = new Int_t  [npart];
 716     fPtT       = new Float_t[npart];
 717     fEtaT      = new Float_t[npart];
 718     fPhiT      = new Float_t[npart];
 719     fPdgT      = new Int_t[npart];
 720
 721     fNt   = npart;
 722     fNtS  = 0;
 723     Float_t chTmp=0.0; // charge of current particle
 724     //    Int_t idGeant;
 725
 726     // this is for Pythia ??
 727     for (Int_t part = 0; part < npart; part++) {
 728         TParticle *MPart = gAlice->Particle(part);
 729         Int_t mpart   = MPart->GetPdgCode();
 730         Int_t child1  = MPart->GetFirstDaughter();
 731         Float_t pT    = MPart->Pt();
 732         Float_t p     = MPart->P();
 733         Float_t phi   = MPart->Phi();
 734         Float_t eta   = -100.;
 735         if(pT > 0.001) eta   = MPart->Eta();
 736         Float_t theta = MPart->Theta();
 737         if  (fDebug>=2) {
 738            printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
 739            part, mpart, child1, MPart->GetLastDaughter(), MPart->GetStatusCode());
 740         }
 741
 742         if (fDebug >= 2) {
 743             if (part == 6 || part == 7)
 744             {
 745                 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
 746                        part-5, pT, eta, phi);
 747             }
 748
 749 //          if (mpart == 21)
 750
 751 //              printf("\n Simulated Jet (pt, eta, phi): %d %d %f %f %f",
 752 //                     part, mpart, pT, eta, phi);
 753         }
 754
 755         fTrackList[part] = 0;
 756         fPtT[part]       = pT; // must be change after correction for resolution !!!
 757         fEtaT[part]      = eta;
 758         fPhiT[part]      = phi;
 759         fPdgT[part]      = mpart;
 760
 761
 762         if (part < 2) continue;
 763
 764         // move to fLego->Fill because hadron correction must apply
 765         // if particle hit to EMCAL - 11-feb-2002
 766         //      if (pT == 0 || pT < fPtCut) continue;
 767         TParticlePDG* pdgP = 0;
 768 // charged or neutral
 769         pdgP  = MPart->GetPDG();
 770         chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
 771
 772         if (ich == 0) {
 773             if (chTmp == 0) {
 774                 if (!fK0N) {
 775                     continue;
 776                 } else {
 777                     if (mpart != kNeutron    &&
 778                         mpart != kNeutronBar &&
 779                         mpart != kK0Long) continue;
 780                 }
 781             }
 782         }
 783
 784 // skip partons
 785         if (TMath::Abs(mpart) <= 6         ||
 786             mpart == 21                    ||
 787             mpart == 92) continue;
 788
 789         if (TMath::Abs(eta)<=0.9) fNChTpc++;
 790 // final state only
 791         if (child1 >= 0 && child1 < npart) continue;
 792 // acceptance cut
 793         if (eta > fEtaMax || eta < fEtaMin)    continue;
 794         if (phi > fPhiMax || phi < fPhiMin )   continue;
 795 //
 796         if (fDebug >= 3)
 797         printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
 798         part, mpart, child1, eta, phi, pT, chTmp);
 799 //
 800 //
 801 // Momentum smearing goes here ...
 802 //
 803         fhTrackPt->Fill(pT);
 804         Float_t pw;
 805         if (fSmear && TMath::Abs(chTmp)) {
 806             pw = AliEMCALFast::SmearMomentum(1,p);
 807         // p changed - take into account when calculate pT,
 808         // pz and so on ..  ?
 809             pT = (pw/p) * pT;
 810             if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
 811             p  = pw;
 812         }
 813 //
 814 // Tracking Efficiency and TPC acceptance goes here ...
 815         Float_t eff;
 816         if (fEffic && TMath::Abs(chTmp)) {
 817           //        eff =  AliEMCALFast::Efficiency(1,p);
 818             eff = 0.95; // for testing 25-feb-2002
 819             if(fhEff) fhEff->Fill(p, eff);
 820             if (AliEMCALFast::RandomReject(eff)) {
 821               if(fDebug >= 5) printf(" reject due to unefficiency ");
 822               continue;
 823             }
 824         }
 825 //
 826 // Correction of Hadronic Energy goes here ...
 827 //
 828 //
 829 // phi propagation for hadronic correction
 830
 831         Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
 832         Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
 833         if(TMath::Abs(chTmp)) {
 834         // hadr. correction only for charge particle
 835           dphi  = PropagatePhi(pT, chTmp, curls);
 836           phiHC = phi + dphi;
 837           if (fDebug >= 6) {
 838              printf("\n Delta phi %f pT %f ", dphi, pT);
 839              if (curls) printf("\n !! Track is curling");
 840           }
 841           if(!curls) fhTrackPtBcut->Fill(pT);
 842
 843           if (fHCorrection && !curls) {
 844               if (!fHadronCorrector)
 845                  Fatal("AliEMCALJetFinder",
 846                     "Hadronic energy correction required but not defined !");
 847
 848               dpH    = fHadronCorrector->GetEnergy(p, eta, 7);
 849               eTdpH  = dpH*TMath::Sin(theta);
 850
 851               if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n",
 852               phi, phiHC, -eTdpH); // correction is negative
 853               fLego->Fill(eta, phiHC, -eTdpH);
 854               fhLegoHadrCorr->Fill(eta, phiHC, eTdpH);
 855           }
 856         }
 857 //
 858 //  More to do ? Just think about it !
 859 //
 860         if (phi > fPhiMax || phi < fPhiMin )   continue;
 861
 862         if(TMath::Abs(chTmp) ) { // charge particle
 863           if (pT > fPtCut && !curls) {
 864              if (fDebug >= 8) printf("Charge :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
 865                                       eta , phi, pT, fNtS);
 866              fLego->Fill(eta, phi, pT);
 867              fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
 868              fTrackList[part] = 1;
 869              fNtS++;
 870           }
 871         } else if(ich==0 && fK0N) {
 872           // case of n, nbar and K0L
 873              if (fDebug >= 9) printf("Neutral :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
 874                                       eta , phi, pT, fNtS);
 875             fLego->Fill(eta, phi, pT);
 876             fTrackList[part] = 1;
 877             fNtS++;
 878         }
 879
 880     } // primary loop
 881     DumpLego();
 882 }
 883
 884 void AliEMCALJetFinder::FillFromHits(Int_t flag)
 885 {
 886 //
 887 // Fill Cells with hit information
 888 //
 889 //
 890     if (fDebug >= 2)
 891     printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
 892
 893     ResetMap();
 894
 895     if (!fLego) BookLego();
 896 //  Reset eta-phi maps if needed
 897     if (flag == 0) { // default behavior
 898       fLego->Reset();
 899       fhLegoTracks->Reset();
 900       fhLegoEMCAL->Reset();
 901       fhLegoHadrCorr->Reset();
 902     }
 903 //  Initialize from background event if available
 904 //
 905 // Access hit information
 906     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 907     TTree *treeH = gAlice->TreeH();
 908     Int_t ntracks = (Int_t) treeH->GetEntries();
 909 //
 910 //   Loop over tracks
 911 //
 912     Int_t nbytes = 0;
 913     Double_t etH = 0.0;
 914
 915     for (Int_t track=0; track<ntracks;track++) {
 916         gAlice->ResetHits();
 917         nbytes += treeH->GetEvent(track);
 918 //
 919 //  Loop over hits
 920 //
 921         for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
 922             mHit;
 923             mHit=(AliEMCALHit*) pEMCAL->NextHit())
 924         {
 925             Float_t x      =    mHit->X();         // x-pos of hit
 926             Float_t y      =    mHit->Y();         // y-pos
 927             Float_t z      =    mHit->Z();         // z-pos
 928             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
 929
 930             Float_t r      =    TMath::Sqrt(x*x+y*y);
 931             Float_t theta  =    TMath::ATan2(r,z);
 932             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
 933             Float_t phi    =    TMath::ATan2(y,x);
 934
 935             if (fDebug >= 11) printf("\n Hit %f %f %f %f", x, y, z, eloss);
 936 //          printf("\n Hit %f %f %f %f", x, y, z, eloss);
 937
 938             etH = fSamplingF*eloss*TMath::Sin(theta);
 939             fLego->Fill(eta, phi, etH);
 940             //      fhLegoEMCAL->Fill(eta, phi, etH);
 941         } // Hit Loop
 942     } // Track Loop
 943     // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
 944     for(Int_t i=0; i<fLego->GetSize(); i++) (*fhLegoEMCAL)[i] = (*fLego)[i];
 945     //    DumpLego(); ??
 946 }
 947
 948 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
 949 {
 950 //
 951 // Fill Cells with digit information
 952 //
 953 //
 954     ResetMap();
 955
 956     if (!fLego) BookLego();
 957     if (flag == 0) fLego->Reset();
 958     Int_t nbytes;
 959
 960
 961 //
 962 //  Connect digits
 963 //
 964     TClonesArray* digs      = new TClonesArray("AliEMCALDigit", 10000);
 965     TTree *treeD = gAlice->TreeD();
 966     TBranchElement* branchDg = (TBranchElement*)
 967         treeD->GetBranch("EMCAL");
 968
 969     if (!branchDg) Fatal("AliEMCALJetFinder",
 970                          "Reading digits requested but no digits in file !");
 971
 972     branchDg->SetAddress(&digs);
 973     Int_t nent = (Int_t) branchDg->GetEntries();
 974 //
 975 //  Connect digitizer
 976 //
 977     AliEMCALDigitizer* digr = new AliEMCALDigitizer();
 978     TBranchElement* branchDr = (TBranchElement*)
 979         treeD->GetBranch("AliEMCALDigitizer");
 980     branchDr->SetAddress(&digr);
 981 //
 982 //
 983     nbytes += branchDg->GetEntry(0);
 984     nbytes += branchDr->GetEntry(0);
 985 //
 986 //  Get digitizer parameters
 987     Float_t preADCped = digr->GetPREpedestal();
 988     Float_t preADCcha = digr->GetPREchannel();
 989     Float_t ecADCped  = digr->GetECpedestal();
 990     Float_t ecADCcha  = digr->GetECchannel();
 991     Float_t hcADCped  = digr->GetHCpedestal();
 992     Float_t hcADCcha  = digr->GetHCchannel();
 993
 994     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 995     AliEMCALGeometry* geom =
 996         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 997
 998     if (fDebug) {
 999         Int_t   ndig = digs->GetEntries();
1000         Info("FillFromDigits","Number of Digits: %d %d\n Parameters: PRE : %f %f EC: %f %f HC: %f %f\n Geometry: %d %d",
1001              ndig, nent, preADCped, preADCcha, ecADCped, ecADCcha, hcADCped, hcADCcha, geom->GetNEta(), geom->GetNPhi());
1002     }
1003
1004 //
1005 //  Loop over digits
1006     AliEMCALDigit* sdg;
1007     TIter next(digs);
1008     while ((sdg = (AliEMCALDigit*)(next())))
1009     {
1010         Double_t pedestal = 0.;
1011         Double_t channel  = 0.;
1012         if (geom->IsInPRE(sdg->GetId())) {
1013           pedestal = preADCped;
1014           channel  = preADCcha;
1015         }
1016         else if (geom->IsInECAL(sdg->GetId()))  {
1017           pedestal = ecADCped;
1018           channel  = ecADCcha;
1019         }
1020         else if (geom->IsInHCAL(sdg->GetId()))  {
1021           pedestal = hcADCped;
1022           channel  = hcADCcha;
1023         }
1024         else
1025           Fatal("FillFromDigits", "unexpected digit is number!") ;
1026
1027         Float_t eta = sdg->GetEta();
1028         Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
1029         Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
1030
1031         if (fDebug > 1)
1032           Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
1033                            eta, phi, amp, sdg->GetAmp());
1034
1035         fLego->Fill(eta, phi, fSamplingF*amp);
1036     } // digit loop
1037 //
1038 //  Dump lego hist
1039     DumpLego();
1040 }
1041
1042
1043 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
1044 {
1045 //
1046 // Fill Cells with hit information
1047 //
1048 //
1049     ResetMap();
1050
1051     if (!fLego) BookLego();
1052 // Reset
1053     if (flag == 0) fLego->Reset();
1054
1055 // Flag charged tracks passing through TPC which
1056 // are associated to EMCAL Hits
1057     BuildTrackFlagTable();
1058
1059 //
1060 // Access particle information
1061     TTree *treeK = gAlice->TreeK();
1062     Int_t ntracks = (Int_t) treeK->GetEntries();
1063
1064     if (fPtT)       delete[] fPtT;
1065     if (fEtaT)      delete[] fEtaT;
1066     if (fPhiT)      delete[] fPhiT;
1067     if (fPdgT)      delete[] fPdgT;
1068
1069     fPtT       = new Float_t[ntracks];
1070     fEtaT      = new Float_t[ntracks];
1071     fPhiT      = new Float_t[ntracks];
1072     fPdgT      = new Int_t[ntracks];
1073
1074     fNt   = ntracks;
1075     fNtS  = 0;
1076
1077     for (Int_t track = 0; track < ntracks; track++) {
1078           TParticle *MPart = gAlice->Particle(track);
1079           Float_t pT    = MPart->Pt();
1080           Float_t phi   = MPart->Phi();
1081           Float_t eta   = MPart->Eta();
1082
1083         if(fTrackList[track]) {
1084           fPtT[track]       = pT;
1085           fEtaT[track]      = eta;
1086           fPhiT[track]      = phi;
1087           fPdgT[track]      = MPart->GetPdgCode();
1088
1089           if (track < 2) continue;      //Colliding particles?
1090           if (pT == 0 || pT < fPtCut) continue;
1091           fNtS++;
1092           fLego->Fill(eta, phi, pT);
1093         }
1094       } // track loop
1095     DumpLego();
1096 }
1097
1098 void AliEMCALJetFinder::FillFromParticles()
1099 {
1100 // 26-feb-2002 PAI - for checking all chain
1101 // Work on particles level; accept all particle (not neutrino )
1102
1103     Double_t PX=0, PY=0, PZ=0, E=0; // checking conservation law
1104     fNChTpc = 0;
1105
1106     ResetMap();
1107     if (!fLego) BookLego();
1108     fLego->Reset();
1109 //
1110 // Access particles information
1111     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1112     if (fDebug >= 2 || npart<=0) {
1113        printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
1114        if(npart<=0) return;
1115     }
1116     fNt   = npart;
1117     fNtS  = 0;
1118     RearrangeParticlesMemory(npart);
1119
1120 //  Go through the particles
1121     Int_t mpart, child1, child2, geantPdg;
1122     Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
1123     TParticle *MPart=0;
1124     for (Int_t part = 0; part < npart; part++) {
1125
1126         fTrackList[part] = 0;
1127
1128         MPart   = gAlice->Particle(part);
1129         mpart   = MPart->GetPdgCode();
1130         child1  = MPart->GetFirstDaughter();
1131         child2  = MPart->GetLastDaughter();
1132         pT      = MPart->Pt();
1133         phi     = MPart->Phi();
1134         eta     = MPart->Eta();
1135
1136         px      = MPart->Px();
1137         py      = MPart->Py();
1138         pz      = MPart->Pz();
1139         e       = MPart->Energy();
1140
1141 // see pyedit in Pythia's text
1142         geantPdg = mpart;
1143         if (IsThisPartonsOrDiQuark(mpart)) continue;
1144         printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n",
1145         part, mpart, geantPdg, px, py, pz, e, child1, child2);
1146
1147 //  exclude partons (21 - gluon, 92 - string)
1148
1149
1150 // exclude neutrinous also ??
1151         if (fDebug >= 11 && pT>0.01)
1152         printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
1153         part, mpart, eta, phi, pT);
1154
1155         fPtT[part]       = pT;
1156         fEtaT[part]      = eta;
1157         fPhiT[part]      = phi;
1158         fPdgT[part]      = mpart;
1159         fNtS++;
1160
1161 // final state only
1162         if (child1 >= 0 && child1 < npart) continue;
1163
1164         //        printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n",
1165         //        part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
1166         PX += px;
1167         PY += py;
1168         PZ += pz;
1169         E  += e;
1170
1171
1172         if (TMath::Abs(eta) <= 0.9) fNChTpc++;
1173 // acceptance cut
1174         if (eta > fEtaMax || eta < fEtaMin)    continue;
1175         if (phi > fPhiMax || phi < fPhiMin )   continue;
1176 //
1177         if(fK0N==0 ) { // exclude neutral hadrons
1178           if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue;
1179         }
1180         fTrackList[part] = 1;
1181         fLego->Fill(eta, phi, pT);
1182
1183     } // primary loop
1184     printf("\n                PX %8.2f  PY %8.2f  PZ %8.2f  E %8.2f \n",
1185     PX, PY, PZ, E);
1186     DumpLego();
1187     if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
1188 }
1189
1190 void AliEMCALJetFinder::FillFromPartons()
1191 {
1192 // function under construction - 13-feb-2002 PAI
1193
1194     if (fDebug >= 2)
1195     printf("\n AliEMCALJetFinder::FillFromPartons()\n");
1196     //
1197
1198     ResetMap();
1199     if (!fLego) BookLego();
1200     fLego->Reset();
1201 //
1202 // Access particle information
1203     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1204     if (fDebug >= 2 || npart<=0)
1205     printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
1206     fNt   = 0; // for FindTracksInJetCone
1207     fNtS  = 0;
1208
1209 //  Go through the partons
1210     Int_t statusCode=0;
1211     for (Int_t part = 8; part < npart; part++) {
1212         TParticle *MPart = gAlice->Particle(part);
1213         Int_t mpart   = MPart->GetPdgCode();
1214         //      Int_t child1  = MPart->GetFirstDaughter();
1215         Float_t pT    = MPart->Pt();
1216         //      Float_t p     = MPart->P();
1217         Float_t phi   = MPart->Phi();
1218         Float_t eta   = MPart->Eta();
1219         //      Float_t theta = MPart->Theta();
1220         statusCode    = MPart->GetStatusCode();
1221
1222 // accept partons (21 - gluon, 92 - string)
1223         if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
1224         if (fDebug > 1 && pT>0.01)
1225         printf("\n part:%5d mpart %5d status  %5d eta %8.2f phi %8.2f pT %8.2f ",
1226         part, mpart, statusCode, eta, phi, pT);
1227         //      if (fDebug >= 3) MPart->Print();
1228 // accept partons before fragmentation - p.57 in Pythia manual
1229 //        if(statusCode != 1) continue;
1230 // acceptance cut
1231         if (eta > fEtaMax || eta < fEtaMin)    continue;
1232         if (phi > fPhiMax || phi < fPhiMin )   continue;
1233 // final state only
1234 //      if (child1 >= 0 && child1 < npart) continue;
1235 //
1236 //
1237         fLego->Fill(eta, phi, pT);
1238
1239     } // primary loop
1240     DumpLego();
1241 }
1242
1243 void AliEMCALJetFinder::BuildTrackFlagTable() {
1244
1245 // Method to generate a lookup table for TreeK particles
1246 // which are linked to hits in the EMCAL
1247 //
1248 // --Author: J.L. Klay
1249 //
1250 // Access hit information
1251     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1252
1253     TTree *TK = gAlice->TreeK();                // Get the number of entries in the kine tree
1254     Int_t nKTrks = (Int_t) TK->GetEntries();    // (Number of particles created somewhere)
1255
1256     if(fTrackList) delete[] fTrackList;         //Make sure we get rid of the old one
1257     fTrackList = new Int_t[nKTrks];             //before generating a new one
1258
1259     for (Int_t i = 0; i < nKTrks; i++) {        //Initialize members to 0
1260         fTrackList[i] = 0;
1261     }
1262
1263     TTree *treeH = gAlice->TreeH();
1264     Int_t ntracks = (Int_t) treeH->GetEntries();
1265 //
1266 //   Loop over tracks
1267 //
1268     Int_t nbytes = 0;
1269
1270     for (Int_t track=0; track<ntracks;track++) {
1271         gAlice->ResetHits();
1272         nbytes += treeH->GetEvent(track);
1273         if (pEMCAL)  {
1274
1275 //
1276 //  Loop over hits
1277 //
1278             for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1279                 mHit;
1280                 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1281             {
1282                 Int_t   iTrk   =    mHit->Track();       // track number
1283                 Int_t   idprim =    mHit->GetPrimary();  // primary particle
1284
1285                 //Determine the origin point of this particle - it made a hit in the EMCAL
1286                 TParticle         *trkPart = gAlice->Particle(iTrk);
1287                 TParticlePDG  *trkPdg  = trkPart->GetPDG();
1288                 Int_t      trkCode = trkPart->GetPdgCode();
1289                 Double_t           trkChg;
1290                 if (trkCode < 10000) {          //Big Ions cause problems for
1291                     trkChg  = trkPdg->Charge(); //this function.  Since they aren't
1292                 } else {                                //likely to make it very far, set
1293                     trkChg  = 0.0;                      //their charge to 0 for the Flag test
1294                 }
1295                 Float_t            vxTrk   = trkPart->Vx();
1296                 Float_t            vyTrk   = trkPart->Vy();
1297                 Float_t            vrTrk   = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
1298                 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1299
1300                 //Loop through the ancestry of the EMCAL entrance particles
1301                 Int_t ancestor = trkPart->GetFirstMother();  //Get track's Mother
1302                 while (ancestor != -1) {
1303                     TParticle     *ancPart = gAlice->Particle(ancestor);  //get particle info on ancestor
1304                     TParticlePDG  *ancPdg  = ancPart->GetPDG();
1305                     Int_t             ancCode = ancPart->GetPdgCode();
1306                     Double_t       ancChg;
1307                     if (ancCode < 10000) {
1308                         ancChg  = ancPdg->Charge();
1309                     } else {
1310                         ancChg  = 0.0;
1311                     }
1312                     Float_t           vxAnc   = ancPart->Vx();
1313                     Float_t           vyAnc   = ancPart->Vy();
1314                     Float_t           vrAnc   = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1315                     fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1316                     ancestor = ancPart->GetFirstMother();           //Get the next ancestor
1317                 }
1318
1319                 //Determine the origin point of the primary particle associated with the hit
1320                 TParticle     *primPart = gAlice->Particle(idprim);
1321                 TParticlePDG  *primPdg  = primPart->GetPDG();
1322                 Int_t      primCode = primPart->GetPdgCode();
1323                 Double_t           primChg;
1324                 if (primCode < 10000) {
1325                     primChg  = primPdg->Charge();
1326                 } else {
1327                     primChg  = 0.0;
1328                 }
1329                 Float_t    vxPrim = primPart->Vx();
1330                 Float_t    vyPrim = primPart->Vy();
1331                 Float_t    vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1332                 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1333             } // Hit Loop
1334         } //if (pEMCAL)
1335     } // Track Loop
1336 }
1337
1338 Int_t AliEMCALJetFinder
1339 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1340
1341     Int_t flag    = 0;
1342     Int_t parton  = 0;
1343     Int_t neutral = 0;
1344
1345     if (charge == 0) neutral = 1;
1346
1347     if (TMath::Abs(code) <= 6   ||
1348         code == 21              ||
1349         code == 92) parton = 1;
1350
1351     //It's not a parton, it's charged and it went through the TPC
1352     if (!parton && !neutral && radius <= 84.0) flag = 1;
1353
1354     return flag;
1355 }
1356
1357
1358
1359 void AliEMCALJetFinder::SaveBackgroundEvent()
1360 {
1361 // Saves the eta-phi lego and the tracklist
1362 //
1363     if (fLegoB) {
1364        fLegoB->Reset();
1365        (*fLegoB) = (*fLegoB) + (*fLego);
1366        if(fDebug)
1367        printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
1368        fLegoB->Integral(), fLego->Integral());
1369     }
1370
1371     if (fPtB)        delete[] fPtB;
1372     if (fEtaB)       delete[] fEtaB;
1373     if (fPhiB)       delete[] fPhiB;
1374     if (fPdgB)       delete[] fPdgB;
1375     if (fTrackListB) delete[] fTrackListB;
1376
1377     fPtB          = new Float_t[fNtS];
1378     fEtaB         = new Float_t[fNtS];
1379     fPhiB         = new Float_t[fNtS];
1380     fPdgB         = new Int_t  [fNtS];
1381     fTrackListB   = new Int_t  [fNtS];
1382
1383     fNtB = 0;
1384
1385     for (Int_t i = 0; i < fNt; i++) {
1386         if (!fTrackList[i]) continue;
1387         fPtB [fNtB]       = fPtT [i];
1388         fEtaB[fNtB]       = fEtaT[i];
1389         fPhiB[fNtB]       = fPhiT[i];
1390         fPdgB[fNtB]       = fPdgT[i];
1391         fTrackListB[fNtB] = 1;
1392         fNtB++;
1393     }
1394     fBackground = 1;
1395     printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
1396 }
1397
1398 void AliEMCALJetFinder::InitFromBackground()
1399 {
1400 //
1401 //
1402     if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
1403
1404     if (fLego) {
1405         fLego->Reset();
1406         (*fLego) = (*fLego) + (*fLegoB);
1407         if(fDebug)
1408             printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
1409                    fLego->Integral(), fLegoB->Integral());
1410     } else {
1411         printf(" => fLego undefined \n");
1412     }
1413 }
1414
1415
1416 void AliEMCALJetFinder::FindTracksInJetCone()
1417 {
1418 //
1419 //  Build list of tracks inside jet cone
1420 //
1421 //  Loop over jets
1422     Int_t njet = Njets();
1423     for (Int_t nj = 0; nj < njet; nj++)
1424     {
1425         Float_t etaj = JetEtaW(nj);
1426         Float_t phij = JetPhiW(nj);
1427         Int_t   nT   = 0;           // number of associated tracks
1428
1429 // Loop over particles in current event
1430         for (Int_t part = 0; part < fNt; part++) {
1431             if (!fTrackList[part]) continue;
1432             Float_t phi      = fPhiT[part];
1433             Float_t eta      = fEtaT[part];
1434             Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1435                                            (phij-phi)*(phij-phi));
1436             if (dr < fConeRadius) {
1437                 fTrackList[part] = nj+2;
1438                 nT++;
1439             } // < ConeRadius ?
1440         } // particle loop
1441
1442 // Same for background event if available
1443         Int_t nTB = 0;
1444         if (fBackground) {
1445             for (Int_t part = 0; part < fNtB; part++) {
1446                 Float_t phi      = fPhiB[part];
1447                 Float_t eta      = fEtaB[part];
1448                 Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1449                                                (phij-phi)*(phij-phi));
1450                 fTrackListB[part] = 1;
1451
1452                 if (dr < fConeRadius) {
1453                     fTrackListB[part] = nj+2;
1454                     nTB++;
1455                 } // < ConeRadius ?
1456             } // particle loop
1457         } // Background available ?
1458
1459         Int_t nT0 = nT + nTB;
1460         printf("Total number of tracks %d\n", nT0);
1461
1462         if (nT0 > 1000) nT0 = 1000;
1463
1464         Float_t* ptT  = new Float_t[nT0];
1465         Float_t* etaT = new Float_t[nT0];
1466         Float_t* phiT = new Float_t[nT0];
1467         Int_t*   pdgT = new Int_t[nT0];
1468
1469         Int_t iT = 0;
1470         Int_t j;
1471
1472         for (Int_t part = 0; part < fNt; part++) {
1473             if (fTrackList[part] == nj+2) {
1474                 Int_t index = 0;
1475                 for (j=iT-1; j>=0; j--) {
1476                     if (fPtT[part] > ptT[j]) {
1477                         index = j+1;
1478                         break;
1479                     }
1480                 }
1481                 for (j=iT-1; j>=index; j--) {
1482                     ptT [j+1]  = ptT [j];
1483                     etaT[j+1]  = etaT[j];
1484                     phiT[j+1]  = phiT[j];
1485                     pdgT[j+1]  = pdgT[j];
1486                 }
1487                 ptT [index] = fPtT [part];
1488                 etaT[index] = fEtaT[part];
1489                 phiT[index] = fPhiT[part];
1490                 pdgT[index] = fPdgT[part];
1491                 iT++;
1492             } // particle associated
1493             if (iT > nT0) break;
1494         } // particle loop
1495
1496         if (fBackground) {
1497             for (Int_t part = 0; part < fNtB; part++) {
1498                 if (fTrackListB[part] == nj+2) {
1499                     Int_t index = 0;
1500                     for (j=iT-1; j>=0; j--) {
1501                         if (fPtB[part] > ptT[j]) {
1502                             index = j+1;
1503
1504                             break;
1505                         }
1506                     }
1507                     for (j=iT-1; j>=index; j--) {
1508                         ptT [j+1]  = ptT [j];
1509                         etaT[j+1]  = etaT[j];
1510                         phiT[j+1]  = phiT[j];
1511                         pdgT[j+1]  = pdgT[j];
1512                     }
1513                     ptT [index] = fPtB [part];
1514                     etaT[index] = fEtaB[part];
1515                     phiT[index] = fPhiB[part];
1516                     pdgT[index] = fPdgB[part];
1517                     iT++;
1518                 } // particle associated
1519                 if (iT > nT0) break;
1520             } // particle loop
1521         } // Background available ?
1522
1523         fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
1524         delete[] ptT;
1525         delete[] etaT;
1526         delete[] phiT;
1527         delete[] pdgT;
1528
1529     } // jet loop loop
1530 }
1531
1532 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1533 {
1534 // Propagates phi angle to EMCAL radius
1535 //
1536   static Float_t b = 0.0, rEMCAL = -1.0;
1537   if(rEMCAL<0) {
1538 // Get field in kGS
1539     b =  gAlice->Field()->SolenoidField();
1540 // Get EMCAL radius in cm
1541     rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1542     printf("\nMagnetic field %f rEMCAL %f ", b, rEMCAL);
1543   }
1544     Float_t dPhi = 0.;
1545 //
1546 //
1547 // bending radies
1548 // pt [Gev]
1549 // B  [kG]
1550 //
1551     Float_t rB = 3335.6 * pt / b;  // [cm]  (case of |charge|=1)
1552 //
1553 // check if particle is curling below EMCAL
1554     if (2.*rB < rEMCAL) {
1555         curls = kTRUE;
1556         return dPhi;
1557     }
1558 //
1559 // if not calculate delta phi
1560     Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1561     dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1562     dPhi = -TMath::Sign(dPhi, charge);
1563 //
1564     return dPhi;
1565 }
1566
1567 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
1568 {
1569 // dummy for hbook calls
1570     ;
1571 }
1572
1573 void AliEMCALJetFinder::DrawLego(Char_t *opt)
1574 {fLego->Draw(opt);}
1575
1576 void AliEMCALJetFinder::DrawLegoEMCAL(Char_t *opt)
1577 {fhLegoEMCAL->Draw(opt);}
1578
1579 void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1580 {
1581   static TPaveText *varLabel=0;
1582   if(!fC1) {
1583     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1584   }
1585   fC1->Clear();
1586   fC1->Divide(2,2);
1587   fC1->cd(1);
1588   gPad->SetLogy(1);
1589   fhCellEt->Draw();
1590
1591   fC1->cd(2);
1592   gPad->SetLogy(1);
1593   fhCellEMCALEt->Draw();
1594
1595   fC1->cd(3);
1596   gPad->SetLogy(1);
1597   fhTrackPt->Draw();
1598   fhTrackPtBcut->SetLineColor(2);
1599   fhTrackPtBcut->Draw("same");
1600
1601   fC1->cd(4);
1602   if(!varLabel) {
1603     PrintParameters(1);
1604     varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1605     varLabel->SetTextAlign(12);
1606     varLabel->SetFillColor(19); // see TAttFill
1607     TString tmp(GetTitle());
1608     varLabel->ReadFile(GetFileNameForParameters());
1609   }
1610   varLabel->Draw();
1611   fC1->Update();
1612   if(mode) { // for saving picture to the file
1613     TString stmp(GetFileNameForParameters());
1614     stmp.ReplaceAll("_Par.txt",".ps");
1615     fC1->Print(stmp.Data());
1616   }
1617 }
1618
1619 void AliEMCALJetFinder::PrintParameters(Int_t mode)
1620 {
1621   FILE *file=0;
1622   if(mode==0) file = stdout; // output to terminal
1623   else {
1624     file = fopen(GetFileNameForParameters(),"w");
1625     if(file==0) file = stdout;
1626   }
1627   fprintf(file,"====   Filling lego   ==== \n");
1628   fprintf(file,"Smearing          %6i  ", fSmear);
1629   fprintf(file,"Efficiency        %6i\n", fEffic);
1630   fprintf(file,"Hadr.Correct.     %6i  ", fHCorrection);
1631   fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1632   fprintf(file,"====  Jet finding     ==== \n");
1633   fprintf(file,"Cone radius       %6.2f  ", fConeRadius);
1634   fprintf(file,"Seed E_{T}           %6.1f\n", fEtSeed);
1635   fprintf(file,"Min E_{T} of cell    %6.1f  ", fMinCellEt);
1636   fprintf(file,"Min E_{T} of jet     %6.1f\n", fMinJetEt);
1637   if(fMode) {
1638     fprintf(file,"====  Bg subtraction     ==== \n");
1639     fprintf(file,"BG subtraction  %6i  ", fMode);
1640     fprintf(file,"Min cone move   %6.2f\n", fMinMove);
1641     fprintf(file,"Max cone move   %6.2f  ", fMaxMove);
1642     fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1643   } else
1644   fprintf(file,"==== No Bg subtraction     ==== \n");
1645   if(file != stdout) fclose(file);
1646 }
1647
1648 void AliEMCALJetFinder::DrawLegos()
1649 {
1650   if(!fC1) {
1651     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1652   }
1653   fC1->Clear();
1654   fC1->Divide(2,2);
1655   gStyle->SetOptStat(111111);
1656
1657   Int_t nent1, nent2, nent3, nent4;
1658   Double_t int1, int2, int3, int4;
1659   nent1 = (Int_t)fLego->GetEntries();
1660   int1  = fLego->Integral();
1661   fC1->cd(1);
1662   if(int1) fLego->Draw("lego");
1663
1664   nent2 = (Int_t)fhLegoTracks->GetEntries();
1665   int2  = fhLegoTracks->Integral();
1666   fC1->cd(2);
1667   if(int2) fhLegoTracks->Draw("lego");
1668
1669   nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1670   int3  = fhLegoEMCAL->Integral();
1671   fC1->cd(3);
1672   if(int3) fhLegoEMCAL->Draw("lego");
1673
1674   nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1675   int4  = fhLegoHadrCorr->Integral();
1676   fC1->cd(4);
1677   if(int4) fhLegoHadrCorr->Draw("lego");
1678
1679   // just for checking
1680   printf(" Integrals \n");
1681   printf("lego   %10.3f \ntrack  %10.3f \nhits   %10.3f \nHCorr   %10.3f\n--      %10.3f(must be 0)\n",
1682   int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1683 }
1684
1685 const Char_t* AliEMCALJetFinder::GetFileNameForParameters(Char_t* dir)
1686 {
1687   static TString tmp;
1688   if(strlen(dir)) tmp = dir;
1689   tmp += GetTitle();
1690   tmp += "_Par.txt";
1691   return tmp.Data();
1692 }
1693
1694 void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1695 { // See FillFromTracks() - npart must be positive
1696     if (fTrackList) delete[] fTrackList;
1697     if (fPtT)       delete[] fPtT;
1698     if (fEtaT)      delete[] fEtaT;
1699     if (fPhiT)      delete[] fPhiT;
1700     if (fPdgT)      delete[] fPdgT;
1701
1702     if(npart>0) {
1703        fTrackList = new Int_t  [npart];
1704        fPtT       = new Float_t[npart];
1705        fEtaT      = new Float_t[npart];
1706        fPhiT      = new Float_t[npart];
1707        fPdgT      = new Int_t[npart];
1708     } else {
1709        printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1710     }
1711 }
1712
1713 Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1714 {
1715   Int_t absPdg = TMath::Abs(pdg);
1716   if(absPdg<=6) return kTRUE; // quarks
1717   if(pdg == 21) return kTRUE; // gluon
1718   if(pdg == 92) return kTRUE; // string
1719
1720   // see p.51 of Pythia Manual
1721   // Not include diquarks with c and b quark - 4-mar-2002
1722   //                 ud_0,sd_0,su_0; dd_1,ud_1,uu_1;  sd_1,su_1,ss_1
1723   static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203,  3103,3203,3303};
1724   for(Int_t i=0; i<9; i++) if(absPdg == diquark[i])  return kTRUE; // diquarks
1725
1726   return kFALSE;
1727 }
1728
1729 void AliEMCALJetFinder::FindChargedJet()
1730 {
1731 //
1732 // Find jet from charged particle information only
1733 //
1734
1735 //
1736 //  Look for seeds
1737 //
1738     Int_t njets = 0;
1739     Int_t part  = 0;
1740     Int_t nseed = 0;
1741
1742 //
1743 //
1744     ResetJets();
1745
1746 //
1747     for (part = 0; part < fNt; part++) {
1748         if (!fTrackList[part]) continue;
1749         if (fPtT[part] > fEtSeed) nseed++;
1750     }
1751     printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
1752     Int_t* iSeeds = new Int_t[nseed];
1753     nseed = 0;
1754
1755     for (part = 0; part < fNt; part++) {
1756         if (!fTrackList[part]) continue;
1757         if (fPtT[part] > fEtSeed) iSeeds[nseed++] =  part;
1758     }
1759
1760 //
1761 // Loop over seeds
1762 //
1763     Int_t seed = 0;
1764     Float_t pt;
1765
1766     while(1){
1767 //
1768 // Find seed with highest pt
1769 //
1770         Float_t ptmax = -1.;
1771         Int_t   index = -1;
1772         Int_t   jndex = -1;
1773         for (seed = 0; seed < nseed; seed++) {
1774             if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
1775                 ptmax = pt;
1776                 index = seed;
1777             } // ptmax ?
1778         } // seeds
1779         if (ptmax < 0.) break;
1780         jndex = iSeeds[index];
1781 //
1782 // Remove from the list
1783         iSeeds[index] = -1;
1784         printf("\n Next Seed %d %f", jndex, ptmax);
1785 //
1786 // Find tracks in cone around seed
1787 //
1788         Float_t phiSeed = fPhiT[jndex];
1789         Float_t etaSeed = fEtaT[jndex];
1790         Float_t eT = 0.;
1791         Float_t pxJ = 0.;
1792         Float_t pyJ = 0.;
1793         Float_t pzJ = 0.;
1794
1795         for (part = 0; part < fNt; part++) {
1796             if (!fTrackList[part]) continue;
1797             Float_t deta = fEtaT[part] - etaSeed;
1798             Float_t dphi = fPhiT[part] - phiSeed;
1799             Float_t dR   = TMath::Sqrt(deta * deta + dphi * dphi);
1800             if (dR < fConeRadius) {
1801                 eT += fPtT[part];
1802                 Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
1803                 Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
1804                 Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
1805                 Float_t pz = fPtT[part] / TMath::Tan(theta);
1806                 pxJ += px;
1807                 pyJ += py;
1808                 pzJ += pz;
1809                 //
1810                 // if seed, remove it
1811                 //
1812                 for (seed = 0; seed < nseed; seed++) {
1813                     if (part == iSeeds[seed]) iSeeds[seed] = -1;
1814                 } // seed ?
1815             } // < cone radius
1816         } // particle loop
1817
1818 //
1819 //      Estimate of jet direction
1820         Float_t phiJ   = TMath::ATan2(pyJ, pxJ);
1821         Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
1822         Float_t etaJ   = TMath::Log(TMath::Tan(thetaJ / 2.));
1823 //      Float_t ptJ    = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
1824
1825 //
1826 //      Sum up all energy
1827 //
1828         Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
1829         Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
1830         Int_t dIphi = Int_t(fConeRadius / fDphi);
1831         Int_t dIeta = Int_t(fConeRadius / fDeta);
1832         Int_t iPhi, iEta;
1833         Float_t sumE = 0;
1834         for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
1835             for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
1836                 if (iPhi < 0 || iEta < 0) continue;
1837                 Float_t dPhi = fPhiMin + iPhi * fDphi;
1838                 Float_t dEta = fEtaMin + iEta * fDeta;
1839                 if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
1840                 sumE += fLego->GetBinContent(iEta, iPhi);
1841             } // eta
1842         } // phi
1843 //
1844 //
1845 //
1846         fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
1847         FindTracksInJetCone();
1848         printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
1849         printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
1850         printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
1851     } // while(1)
1852     EMCALJETS.njet = njets;
1853     if (fWrite) WriteJets();
1854     fEvent++;
1855 }