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