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