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