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