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