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