EMCAL/jetfinder/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 /* $Id$ */
  17
  18 //*-- Authors: Andreas Morsch   (CERN)
  19 //*            J.L. Klay        (LBL)
  20 //*            Aleksei Pavlinov (WSU)
  21 //--
  22 //--
  23 //
  24
  25 #include <stdio.h>
  26 // From root ...
  27 #include <TArrayF.h>
  28 #include <TAxis.h>
  29 #include <TBranchElement.h>
  30 #include <TCanvas.h>
  31 #include <TClonesArray.h>
  32 #include <TFile.h>
  33 #include <TH1.h>
  34 #include <TH2.h>
  35 #include <TList.h>
  36 #include <TPDGCode.h>
  37 #include <TPad.h>
  38 #include <TParticle.h>
  39 #include <TParticlePDG.h>
  40 #include <TPaveText.h>
  41 #include <TPythia6Calls.h>
  42 #include <TROOT.h>
  43 #include <TStyle.h>
  44 #include <TTree.h>
  45 #include <TBrowser.h>
  46
  47 // From AliRoot ...
  48 #include "AliEMCALJetFinder.h"
  49 #include "AliHeader.h"
  50 #include "AliMagF.h"
  51 #include "AliMagFCM.h"
  52 #include "AliRun.h"
  53 #include "AliGenerator.h"
  54 #include "AliRunLoader.h"
  55 #include "AliEMCAL.h"
  56 #include "AliEMCALLoader.h"
  57 #include "AliEMCALDigit.h"
  58 #include "AliEMCALDigitizer.h"
  59 #include "AliEMCALFast.h"
  60 #include "AliEMCALGeometry.h"
  61 #include "AliEMCALHadronCorrection.h"
  62 #include "AliEMCALHit.h"
  63 #include "AliEMCALJetMicroDst.h"
  64
  65 // Interface to FORTRAN
  66 #include "Ecommon.h"
  67 #include "AliMC.h"
  68
  69
  70 ClassImp(AliEMCALJetFinder)
  71
  72 //____________________________________________________________________________
  73 AliEMCALJetFinder::AliEMCALJetFinder()
  74   : fBGFileName(""), fEMCALWeight(0.), fTrackWeight(0.), fRandomBg(kFALSE),
  75     fWrite(kTRUE), fWeightingMethod(kFALSE), fJets(0), fLego(0), fLegoB(0),
  76     fhLegoTracks(0), fhLegoEMCAL(0), fhLegoHadrCorr(0), fhEff(0), fhCellEt(0),
  77     fhCellEMCALEt(0), fhTrackPt(0), fhTrackPtBcut(0), fhChPartMultInTpc(0),
  78     fhSinTheta(0), fC1(0), fHistsList(0), fHadronCorrector(0), fHCorrection(0),
  79     fDebug(0), fBackground(0), fConeRadius(0.3),
  80     fPtCut(0.), fEtSeed(0.), fMinJetEt(0.), fMinCellEt(0.),
  81     fSamplingF(0.), fSmear(0), fEffic(0), fK0N(0),
  82     fNjets(0), fDeta(0.), fDphi(0.), fNcell(0), fNtot(0),
  83     fNbinEta(0), fNbinPhi(0), fEtaMin(0.), fEtaMax(0.), fPhiMin(0.),
  84     fPhiMax(0.), fNt(0), fNChTpc(0), fNtS(0), fTrackList(0), fPtT(0),
  85     fEtaT(0), fPhiT(0), fPdgT(0), fNtB(0), fTrackListB(0), fPtB(0),
  86     fEtaB(0),fPhiB(0), fPdgB(0), fMode(0), fMinMove(0.), fMaxMove(0.),
  87     fPrecBg(0.), fError(0), fOutFileName(0), fOutFile(0), fInFile(0), fEvent(0)
  88 {
  89 // Default constructor
  90     fJets             = 0;
  91     fNjets            = 0;
  92     fLego             = 0;
  93     fLegoB            = 0;
  94
  95     fTrackList        = 0;
  96     fPtT              = 0;
  97     fEtaT             = 0;
  98     fPhiT             = 0;
  99     fPdgT             = 0;
 100
 101     fTrackListB       = 0;
 102     fPtB              = 0;
 103     fEtaB             = 0;
 104     fPhiB             = 0;
 105     fPdgB             = 0;
 106
 107     fHCorrection      = 0;
 108     fHadronCorrector  = 0;
 109
 110     fWrite            = 0;
 111     fOutFileName      = 0;
 112     fOutFile          = 0;
 113     fInFile           = 0;
 114     fEvent            = 0;
 115
 116     fRandomBg         = 0;
 117
 118     SetParametersForBgSubtraction();
 119 }
 120
 121 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
 122   : TTask(name, title),
 123     fBGFileName(""), fEMCALWeight(0.), fTrackWeight(0.), fRandomBg(kFALSE),
 124     fWrite(kTRUE), fWeightingMethod(kFALSE), fJets(0), fLego(0), fLegoB(0),
 125     fhLegoTracks(0), fhLegoEMCAL(0), fhLegoHadrCorr(0), fhEff(0), fhCellEt(0),
 126     fhCellEMCALEt(0), fhTrackPt(0), fhTrackPtBcut(0), fhChPartMultInTpc(0),
 127     fhSinTheta(0), fC1(0), fHistsList(0), fHadronCorrector(0), fHCorrection(0),
 128     fDebug(0), fBackground(0), fConeRadius(0.3),
 129     fPtCut(0.), fEtSeed(0.), fMinJetEt(0.), fMinCellEt(0.),
 130     fSamplingF(0.), fSmear(0), fEffic(0), fK0N(0),
 131     fNjets(0), fDeta(0.), fDphi(0.), fNcell(0), fNtot(0),
 132     fNbinEta(0), fNbinPhi(0), fEtaMin(0.), fEtaMax(0.), fPhiMin(0.),
 133     fPhiMax(0.), fNt(0), fNChTpc(0), fNtS(0), fTrackList(0), fPtT(0),
 134     fEtaT(0), fPhiT(0), fPdgT(0), fNtB(0), fTrackListB(0), fPtB(0),
 135     fEtaB(0),fPhiB(0), fPdgB(0), fMode(0), fMinMove(0.), fMaxMove(0.),
 136     fPrecBg(0.), fError(0), fOutFileName(0), fOutFile(0), fInFile(0), fEvent(0)
 137 {
 138 // Constructor
 139 // Title is used in method GetFileNameForParameters();
 140 //
 141     fJets  = new TClonesArray("AliEMCALJet",10000);
 142     fNjets = 0;
 143     for (Int_t i = 0; i < 30000; i++)
 144     {
 145         fEtCell[i]  = 0.;
 146         fEtaCell[i] = 0.;
 147         fPhiCell[i] = 0.;
 148     }
 149     fLego       = 0;
 150     fLegoB      = 0;
 151
 152     fTrackList  = 0;
 153     fPtT        = 0;
 154     fEtaT       = 0;
 155     fPhiT       = 0;
 156     fPdgT       = 0;
 157
 158     fTrackListB       = 0;
 159     fPtB        = 0;
 160     fEtaB       = 0;
 161     fPhiB       = 0;
 162     fPdgB       = 0;
 163
 164     fHCorrection      = 0;
 165     fHadronCorrector  = 0;
 166     fBackground       = 0;
 167     fWrite            = 0;
 168     fOutFileName      = 0;
 169     fOutFile          = 0;
 170     fInFile           = 0;
 171     fEvent            = 0;
 172 //
 173     SetPtCut();
 174     SetMomentumSmearing();
 175     SetEfficiencySim();
 176     SetDebug();
 177     SetHadronCorrection();
 178     SetIncludeK0andN();
 179
 180     fRandomBg         = 0;
 181     SetParametersForBgSubtraction();
 182 }
 183
 184
 185 AliEMCALJetFinder::AliEMCALJetFinder(const AliEMCALJetFinder& jf)
 186   : TTask(jf.GetName(), jf.GetTitle()),
 187     fBGFileName(jf.fBGFileName), fEMCALWeight(jf.fEMCALWeight), fTrackWeight(jf.fTrackWeight),
 188     fRandomBg(jf.fRandomBg),fWrite(jf.fWrite), fWeightingMethod(jf.fWeightingMethod), fJets(jf.fJets),
 189     fLego(jf.fLego), fLegoB(jf.fLegoB), fhLegoTracks(jf.fhLegoTracks), fhLegoEMCAL(jf.fhLegoEMCAL),
 190     fhLegoHadrCorr(jf.fhLegoHadrCorr), fhEff(jf.fhEff), fhCellEt(jf.fhCellEt),fhCellEMCALEt(jf.fhCellEMCALEt),
 191     fhTrackPt(jf.fhTrackPt), fhTrackPtBcut(jf.fhTrackPtBcut), fhChPartMultInTpc(jf.fhChPartMultInTpc),
 192     fhSinTheta(jf.fhSinTheta), fC1(jf.fC1), fHistsList(jf.fHistsList), fHadronCorrector(jf.fHadronCorrector),
 193     fHCorrection(jf.fHCorrection),fDebug(jf.fDebug), fBackground(jf.fBackground), fConeRadius(jf.fConeRadius),
 194     fPtCut(jf.fPtCut), fEtSeed(jf.fEtSeed), fMinJetEt(jf.fMinJetEt), fMinCellEt(jf.fMinCellEt),
 195     fSamplingF(jf.fSamplingF), fSmear(jf.fSmear), fEffic(jf.fEffic), fK0N(jf.fK0N), fNjets(jf.fNjets),
 196     fDeta(jf.fDeta), fDphi(jf.fDphi), fNcell(jf.fNcell), fNtot(jf.fNtot), fNbinEta(jf.fNbinEta),
 197     fNbinPhi(jf.fNbinPhi), fEtaMin(jf.fEtaMin), fEtaMax(jf.fEtaMax), fPhiMin(jf.fPhiMin),
 198     fPhiMax(jf.fPhiMax), fNt(jf.fNt), fNChTpc(jf.fNChTpc), fNtS(jf.fNtS), fTrackList(jf.fTrackList),
 199     fPtT(jf.fPtT), fEtaT(jf.fEtaT), fPhiT(jf.fPhiT), fPdgT(jf.fPdgT), fNtB(jf.fNtB), fTrackListB(jf.fTrackListB),
 200     fPtB(jf.fPtB), fEtaB(jf.fEtaB),fPhiB(jf.fPhiB), fPdgB(jf.fPdgB), fMode(jf.fMode), fMinMove(jf.fMinMove),
 201     fMaxMove(jf.fMaxMove), fPrecBg(jf.fPrecBg), fError(jf.fError), fOutFileName(jf.fOutFileName),
 202     fOutFile(jf.fOutFile), fInFile(jf.fInFile), fEvent(jf.fEvent)
 203 {
 204 // Copy Constructor
 205 }
 206
 207 void AliEMCALJetFinder::SetParametersForBgSubtraction
 208 (Int_t mode, Float_t minMove, Float_t maxMove, Float_t precBg)
 209 {
 210 // see file /home/pavlinov/cosmic/pythia/jetFinderParamData.inc
 211 // at WSU Linux cluster - 11-feb-2002
 212 // These parameters must be tuned more carefull !!!
 213   SetMode(mode);
 214   SetMinMove(minMove);
 215   SetMaxMove(maxMove);
 216   SetPrecBg(precBg);
 217 }
 218
 219 //____________________________________________________________________________
 220 AliEMCALJetFinder::~AliEMCALJetFinder()
 221 {
 222 // Destructor
 223 //
 224     if (fJets){
 225         fJets->Delete();
 226         delete fJets;
 227     }
 228     delete fLego;
 229     delete fLegoB;
 230     delete fhLegoTracks;
 231     delete fhLegoEMCAL;
 232     delete fhLegoHadrCorr;
 233     delete fhEff;
 234     delete fhCellEt;
 235     delete fhCellEMCALEt;
 236     delete fhTrackPt;
 237     delete fhTrackPtBcut;
 238     delete fhChPartMultInTpc;
 239
 240     delete[] fTrackList;
 241     delete[] fPtT;
 242     delete[] fEtaT;
 243     delete[] fPhiT;
 244     delete[] fPdgT;
 245
 246     delete[] fTrackListB;
 247     delete[] fPtB;
 248     delete[] fEtaB;
 249     delete[] fPhiB;
 250     delete[] fPdgB;
 251 }
 252
 253 #ifndef WIN32
 254 # define jet_finder_ua1 jet_finder_ua1_
 255 # define hf1 hf1_
 256 # define type_of_call
 257
 258 #else
 259 # define jet_finder_ua1 JET_FINDER_UA1
 260 # define hf1 HF1
 261 # define type_of_call _stdcall
 262 #endif
 263
 264 extern "C" void type_of_call
 265 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
 266                Float_t  etc[30000],  Float_t etac[30000],
 267                Float_t  phic[30000],
 268                Float_t& min_move, Float_t& max_move, Int_t& mode,
 269                Float_t& prec_bg,  Int_t& ierror);
 270
 271 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
 272
 273
 274 void AliEMCALJetFinder::Init()
 275 {
 276 //
 277 // Geometry and I/O initialization
 278 //
 279 //
 280 //
 281 //  Get geometry parameters from EMCAL
 282 //
 283 //
 284 //  Geometry
 285   //AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 286   AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance();
 287
 288 //    SetSamplingFraction(geom->GetSampling());
 289
 290     fNbinEta = geom->GetNZ();
 291     fNbinPhi = geom->GetNPhi();
 292     fPhiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
 293     fPhiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
 294     fEtaMin  = geom->GetArm1EtaMin();
 295     fEtaMax  = geom->GetArm1EtaMax();
 296     fDphi    = (fPhiMax-fPhiMin)/fNbinPhi;
 297     fDeta    = (fEtaMax-fEtaMin)/fNbinEta;
 298     fNtot    = fNbinPhi*fNbinEta;
 299     fWeightingMethod = kFALSE;
 300
 301 //
 302     SetCellSize(fDeta, fDphi);
 303 //
 304 //  I/O
 305     if (fOutFileName) fOutFile = new TFile(fOutFileName, "recreate");
 306 //
 307 //
 308 }
 309
 310 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncelltot, Float_t etc[30000],
 311                              Float_t etac[30000], Float_t phic[30000],
 312                              Float_t minmove, Float_t maxmove, Int_t mode,
 313                              Float_t precbg,  Int_t   ierror)
 314 {
 315 // Wrapper for fortran coded jet finder
 316 // Full argument list
 317     jet_finder_ua1(ncell, ncelltot, etc, etac, phic,
 318                    minmove, maxmove, mode, precbg, ierror);
 319     // Write result to output
 320     if(fWrite) WriteJets();
 321     fEvent++;
 322 }
 323
 324 void AliEMCALJetFinder::Find()
 325 {
 326 // Wrapper for fortran coded jet finder using member data for
 327 // argument list
 328
 329     Float_t minmove = fMinMove;
 330     Float_t maxmove = fMaxMove;
 331     Int_t   mode     = fMode;
 332     Float_t precbg  = fPrecBg;
 333     Int_t   ierror;
 334     ResetJets(); // 4-feb-2002 by PAI
 335
 336     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
 337                    minmove, maxmove, mode, precbg, ierror);
 338     fError = ierror;
 339     // Write result to output
 340     Int_t njet = Njets();
 341
 342     for (Int_t nj=0; nj<njet; nj++)
 343     {
 344         if (fWeightingMethod)
 345         {
 346           fJetT[nj] = new AliEMCALJet(WeightedJetEnergy( JetEtaW(nj),JetPhiW(nj) ),
 347                                       JetPhiW(nj),
 348                                       JetEtaW(nj));
 349
 350         }else{
 351
 352         fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
 353                                     JetPhiW(nj),
 354                                     JetEtaW(nj));
 355         }
 356         fJetT[nj]->SetIsWeightedEnergy(fWeightingMethod);
 357         fJetT[nj]->SetEMCALEnergy( EMCALConeEnergy(JetEtaW(nj),JetPhiW(nj)) );
 358         fJetT[nj]->SetTrackEnergy(TrackConeEnergy( JetEtaW(nj),JetPhiW(nj)  ));
 359
 360     }
 361
 362     FindTracksInJetCone();
 363     if(fWrite) WriteJets();
 364     fEvent++;
 365 }
 366
 367
 368 Float_t AliEMCALJetFinder::EMCALConeEnergy(Float_t eta, Float_t phi)
 369 {
 370 // Calculate the energy in the cone
 371 Float_t newenergy = 0.0;
 372 Float_t bineta,binphi;
 373 TAxis *x = fhLegoEMCAL->GetXaxis();
 374 TAxis *y = fhLegoEMCAL->GetYaxis();
 375 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 376         {
 377    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 378       {
 379          bineta = x->GetBinCenter(i);
 380          binphi = y->GetBinCenter(j);
 381        if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 382       {
 383            newenergy += fhLegoEMCAL->GetBinContent(i,j);
 384       }
 385     }
 386 }
 387 return newenergy;
 388 }
 389
 390 Float_t AliEMCALJetFinder::TrackConeEnergy(Float_t eta, Float_t phi)
 391 {
 392 // Calculate the track energy in the cone
 393 Float_t newenergy = 0.0;
 394 Float_t bineta,binphi;
 395 TAxis *x = fhLegoTracks->GetXaxis();
 396 TAxis *y = fhLegoTracks->GetYaxis();
 397 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 398    {
 399    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 400    {
 401     bineta = x->GetBinCenter(i);
 402     binphi = y->GetBinCenter(j);
 403     if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 404       {
 405           newenergy += fhLegoTracks->GetBinContent(i,j);
 406        }
 407     }
 408 }
 409 return newenergy;
 410 }
 411
 412
 413 Float_t AliEMCALJetFinder::WeightedJetEnergy(Float_t eta, Float_t phi)
 414 {
 415 // Calculate the weighted jet energy
 416
 417 Float_t newenergy = 0.0;
 418 Float_t bineta,binphi;
 419 TAxis *x = fhLegoEMCAL->GetXaxis();
 420 TAxis *y = fhLegoEMCAL->GetYaxis();
 421
 422
 423 for (Int_t i = 0 ; i < fNbinEta  ; i++) // x coord
 424 {
 425    for (Int_t j = 0 ; j < fNbinPhi  ; j++)  // y coord
 426    {
 427       bineta = x->GetBinCenter(i);
 428       binphi = y->GetBinCenter(j);
 429       if ( (bineta-eta)*(bineta-eta) + (binphi-phi)*(binphi-phi) < fConeRadius*fConeRadius)
 430       {
 431           newenergy += (fEMCALWeight)* fhLegoEMCAL->GetBinContent(i,j) + (fTrackWeight)* fhLegoTracks->GetBinContent(i,j);
 432       }
 433     }
 434 }
 435
 436 return newenergy;
 437
 438 }
 439
 440
 441 Int_t AliEMCALJetFinder::Njets() const
 442 {
 443 // Get number of reconstructed jets
 444     return EMCALJETS.njet;
 445 }
 446
 447 Float_t AliEMCALJetFinder::JetEnergy(Int_t i) const
 448 {
 449 // Get reconstructed jet energy
 450     return EMCALJETS.etj[i];
 451 }
 452
 453 Float_t AliEMCALJetFinder::JetPhiL(Int_t i) const
 454 {
 455 // Get reconstructed jet phi from leading particle
 456     return EMCALJETS.phij[0][i];
 457 }
 458
 459 Float_t AliEMCALJetFinder::JetPhiW(Int_t i) const
 460 {
 461 // Get reconstructed jet phi from weighting
 462     return EMCALJETS.phij[1][i];
 463 }
 464
 465 Float_t  AliEMCALJetFinder::JetEtaL(Int_t i) const
 466 {
 467 // Get reconstructed jet eta from leading particles
 468     return EMCALJETS.etaj[0][i];
 469 }
 470
 471
 472 Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)  const
 473 {
 474 // Get reconstructed jet eta from weighting
 475     return EMCALJETS.etaj[1][i];
 476 }
 477
 478 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
 479 {
 480 // Set grid cell size
 481     EMCALCELLGEO.etaCellSize = eta;
 482     EMCALCELLGEO.phiCellSize = phi;
 483 }
 484
 485 void AliEMCALJetFinder::SetConeRadius(Float_t par)
 486 {
 487 // Set jet cone radius
 488     EMCALJETPARAM.coneRad = par;
 489     fConeRadius = par;
 490 }
 491
 492 void AliEMCALJetFinder::SetEtSeed(Float_t par)
 493 {
 494 // Set et cut for seeds
 495     EMCALJETPARAM.etSeed = par;
 496     fEtSeed = par;
 497 }
 498
 499 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
 500 {
 501 // Set minimum jet et
 502     EMCALJETPARAM.ejMin = par;
 503     fMinJetEt = par;
 504 }
 505
 506 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
 507 {
 508 // Set et cut per cell
 509     EMCALJETPARAM.etMin = par;
 510     fMinCellEt = par;
 511 }
 512
 513 void AliEMCALJetFinder::SetPtCut(Float_t par)
 514 {
 515 // Set pT cut on charged tracks
 516     fPtCut = par;
 517 }
 518
 519
 520 void AliEMCALJetFinder::Test()
 521 {
 522 //
 523 // Test the finder call
 524 //
 525     const Int_t knmax = 30000;
 526     Int_t ncell      = 10;
 527     Int_t ncelltot  = 100;
 528
 529     Float_t etc[knmax];
 530     Float_t etac[knmax];
 531     Float_t phic[knmax];
 532     Float_t minmove = 0;
 533     Float_t maxmove = 0;
 534     Int_t   mode = 0;
 535     Float_t precbg = 0;
 536     Int_t   ierror = 0;
 537
 538
 539     Find(ncell, ncelltot, etc, etac, phic,
 540          minmove, maxmove, mode, precbg, ierror);
 541
 542 }
 543
 544 //
 545 //  I/O
 546 //
 547
 548 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
 549 {
 550     //
 551     // Add a jet
 552     //
 553     TClonesArray &lrawcl = *fJets;
 554     new(lrawcl[fNjets++]) AliEMCALJet(jet);
 555 }
 556
 557 void AliEMCALJetFinder::ResetJets()
 558 {
 559     //
 560     // Reset Jet List
 561     //
 562     fJets->Clear();
 563     fNjets = 0;
 564 }
 565
 566 void AliEMCALJetFinder::WriteJets()
 567 {
 568 //
 569 // Add all jets to the list
 570 //
 571     const Int_t kBufferSize = 4000;
 572     const char* file = 0;
 573
 574     Int_t njet = Njets();
 575
 576     for (Int_t nj = 0; nj < njet; nj++)
 577     {
 578         AddJet(*fJetT[nj]);
 579         delete fJetT[nj];
 580     }
 581
 582 // I/O
 583     AliRunLoader *rl = AliRunLoader::GetRunLoader();
 584     AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
 585
 586     if (!fOutFileName) {
 587 //
 588 // output written to input file
 589 //
 590       AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
 591       TTree* pK = rl->TreeK();
 592       file = (pK->GetCurrentFile())->GetName();
 593       TBranch * jetBranch ;
 594         if (fDebug > 1)
 595             printf("Make Branch - TreeR address %p %p\n",(void*)gAlice->TreeR(), (void*)pEMCAL);
 596         //if (fJets && gAlice->TreeR()) {
 597         if (fJets && emcalLoader->TreeR()) {
 598           // pEMCAL->MakeBranchInTree(gAlice->TreeR(),
 599           jetBranch = emcalLoader->TreeR()->Branch("EMCALJets", &fJets, kBufferSize, 0) ;
 600           //pEMCAL->MakeBranchInTree(gime->TreeR(),
 601           //                         "EMCALJets",
 602           //                         &fJets,
 603           //                         kBufferSize,
 604           //                         file);
 605
 606           //Int_t nev = gAlice->GetHeader()->GetEvent();
 607           //gAlice->TreeR()->Fill();
 608           jetBranch->Fill();
 609           //char hname[30];
 610           //sprintf(hname,"TreeR%d", nev);
 611           //gAlice->TreeR()->Write(hname);
 612           //gAlice->TreeR()->Reset();
 613           //gime->WriteRecPoints("OVERWRITE");
 614           emcalLoader->WriteRecPoints("OVERWRITE");
 615         }
 616     } else {
 617 //
 618 // Output written to user specified output file
 619 //
 620       //TTree* pK = gAlice->TreeK();
 621       TTree* pK = rl->TreeK();
 622       fInFile  = pK->GetCurrentFile();
 623
 624       fOutFile->cd();
 625       char hname[30];
 626       sprintf(hname,"TreeR%d", fEvent);
 627       TTree* treeJ = new TTree(hname, "EMCALJets");
 628       treeJ->Branch("EMCALJets", &fJets, kBufferSize);
 629       treeJ->Fill();
 630       treeJ->Write(hname);
 631       fInFile->cd();
 632     }
 633     ResetJets();
 634 }
 635
 636 void AliEMCALJetFinder::BookLego()
 637 {
 638 //
 639 //  Book histo for discretization
 640 //
 641
 642 //
 643 //  Don't add histos to the current directory
 644     if(fDebug) printf("\n AliEMCALJetFinder::BookLego() \n");
 645
 646     //    TH2::AddDirectory(0); // hists wil be put to the list from gROOT
 647     //    TH1::AddDirectory(0);
 648     gROOT->cd();
 649 //
 650 //  Signal map
 651     fLego = new TH2F("legoH","eta-phi",
 652                            fNbinEta, fEtaMin, fEtaMax,
 653                            fNbinPhi, fPhiMin, fPhiMax);
 654 //
 655 //  Background map
 656     fLegoB = new TH2F("legoB","eta-phi for BG event",
 657                            fNbinEta, fEtaMin, fEtaMax,
 658                            fNbinPhi, fPhiMin, fPhiMax);
 659
 660 //  Tracks map
 661     fhLegoTracks = new TH2F("hLegoTracks","eta-phi for Tracks",
 662     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 663 //  EMCAL map
 664     fhLegoEMCAL = new TH2F("hLegoEMCAL","eta-phi for EMCAL",
 665     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 666 //  Hadron correction map
 667     fhLegoHadrCorr = new TH2F("hLegoHadrCorr","eta-phi for Hadron. Correction",
 668     fNbinEta, fEtaMin, fEtaMax, fNbinPhi, fPhiMin, fPhiMax);
 669 //  Hists. for tuning jet finder
 670     fhEff = new TH2F("hEff","#epsilon vs momentum ", 100,0.,20., 50,0.5,1.);
 671
 672     TArrayF eTmp(1101);
 673     eTmp[0] = 0.0;
 674     for(Int_t i=1; i<=1000; i++)      eTmp[i] = 0.1*i; // step 100 mev
 675     for(Int_t i=1001; i<=1100; i++)   eTmp[i] = eTmp[1000] + 1.0*(i-1000); // step 1GeV
 676
 677     fhCellEt  = new TH1F("hCellEt","Cell E_{T} from fLego",
 678     eTmp.GetSize()-1, eTmp.GetArray());
 679     fhCellEMCALEt  = new TH1F("hCellEMCALEt","Cell E_{T} for EMCAL itself",
 680     eTmp.GetSize()-1, eTmp.GetArray());
 681     fhTrackPt = new TH1F("hTrackPt","Ch.particles P_{T} ",
 682     eTmp.GetSize()-1, eTmp.GetArray());
 683     fhTrackPtBcut = new TH1F("hTrackPtBcut","Ch.particles P_{T} + magnetic field cut",
 684     eTmp.GetSize()-1, eTmp.GetArray());
 685
 686     fhChPartMultInTpc = new TH1F("hChPartMultInTpc",
 687     "Charge partilcle multiplicity in |%eta|<0.9", 2000, 0, 20000);
 688
 689     fhSinTheta = new TH1F("fhSinTheta","sin(theta)", fNbinEta, fEtaMin, fEtaMax);
 690     TAxis *xax = fhSinTheta->GetXaxis();
 691     for(Int_t i=1; i<=fNbinEta; i++) {
 692       Double_t eta = xax->GetBinCenter(i);
 693       fhSinTheta->Fill(eta, 1./TMath::CosH(eta)); // cosh(eta) = 1./sin(theta)
 694     }
 695
 696     //! first canvas for drawing
 697     fHistsList=AliEMCALJetMicroDst::MoveHistsToList("Hists from AliEMCALJetFinder", kFALSE);
 698 }
 699
 700 void AliEMCALJetFinder::DumpLego()
 701 {
 702 //
 703 // Dump lego histo into array
 704 //
 705     fNcell = 0;
 706     TAxis* xaxis = fLego->GetXaxis();
 707     TAxis* yaxis = fLego->GetYaxis();
 708     //    fhCellEt->Clear();
 709     Float_t e, eH;
 710     for (Int_t i = 1; i <= fNbinEta; i++) {
 711         for (Int_t j = 1; j <= fNbinPhi; j++) {
 712
 713             e = fLego->GetBinContent(i,j);
 714             if (fRandomBg) {
 715                 if (gRandom->Rndm() < 0.5) {
 716                     Float_t ebg = 0.28 * TMath::Abs(gRandom->Gaus(0.,1.));
 717                     e += ebg;
 718                 }
 719             }
 720             if (e > 0.0) e -= fMinCellEt;
 721             if (e < 0.0) e = 0.;
 722             Float_t eta  = xaxis->GetBinCenter(i);
 723             Float_t phi  = yaxis->GetBinCenter(j);
 724             fEtCell[fNcell]  = e;
 725             fEtaCell[fNcell] = eta;
 726             fPhiCell[fNcell] = phi;
 727             fNcell++;
 728             fhCellEt->Fill(e);
 729             if(fhLegoEMCAL) {
 730               eH = fhLegoEMCAL->GetBinContent(i,j);
 731               if(eH > 0.0) fhCellEMCALEt->Fill(eH);
 732             }
 733         } // phi
 734     } // eta
 735 //  today
 736 //    fNcell--;
 737 }
 738
 739 void AliEMCALJetFinder::ResetMap()
 740 {
 741 //
 742 // Reset eta-phi array
 743
 744     for (Int_t i=0; i<30000; i++)
 745     {
 746         fEtCell[i]  = 0.;
 747         fEtaCell[i] = 0.;
 748         fPhiCell[i] = 0.;
 749     }
 750 }
 751
 752
 753 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 754 {
 755 // Which generator
 756 //
 757     const char* name = gAlice->Generator()->GetName();
 758     enum {kPythia, kHijing, kHijingPara};
 759     Int_t genType = 0;
 760
 761     if (!strcmp(name, "Hijing")){
 762         genType = kHijing;
 763     } else if (!strcmp(name, "Pythia")) {
 764         genType = kPythia;
 765     } else if (!strcmp(name, "HIJINGpara") ||!strcmp(name, "HIGINGpara")) {
 766         genType = kHijingPara;
 767     }
 768     if (fDebug>=2)
 769         printf("Fill tracks generated by %s %d\n", name, genType);
 770
 771
 772 //
 773 // Fill Cells with track information
 774 //
 775     if (fDebug >= 2)
 776     printf("\n AliEMCALJetFinder::FillFromTracks(%i,%i) ",flag,ich);
 777     fNChTpc = 0;
 778
 779     ResetMap();
 780
 781     if (!fLego) BookLego();
 782 // Reset
 783     if (flag == 0) fLego->Reset();
 784 //
 785 // Access particle information
 786     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 787     Int_t ntr   = (gAlice->GetHeader())->GetNtrack();
 788     printf(" : #primary particles %i # tracks %i : (before) Sum.Et %f\n",
 789     npart, ntr, fLego->Integral());
 790
 791 // Create track list
 792 //
 793 // 0: not selected
 794 // 1: selected for jet finding
 795 // 2: inside jet
 796 // ....
 797     if (fTrackList) delete[] fTrackList;
 798     if (fPtT)       delete[] fPtT;
 799     if (fEtaT)      delete[] fEtaT;
 800     if (fPhiT)      delete[] fPhiT;
 801     if (fPdgT)      delete[] fPdgT;
 802
 803     fTrackList = new Int_t  [npart];
 804     fPtT       = new Float_t[npart];
 805     fEtaT      = new Float_t[npart];
 806     fPhiT      = new Float_t[npart];
 807     fPdgT      = new Int_t[npart];
 808
 809     fNt   = npart;
 810     fNtS  = 0;
 811     Float_t chTmp=0.0; // charge of current particle
 812     //    Int_t idGeant;
 813
 814     // this is for Pythia ??
 815     for (Int_t part = 0; part < npart; part++) {
 816         TParticle *mPart = gAlice->GetMCApp()->Particle(part);
 817         Int_t mpart   = mPart->GetPdgCode();
 818         Int_t child1  = mPart->GetFirstDaughter();
 819         Float_t pT    = mPart->Pt();
 820         Float_t p     = mPart->P();
 821         Float_t phi   = mPart->Phi();
 822         Float_t eta   = -100.;
 823         if(pT > 0.001) eta   = mPart->Eta();
 824         Float_t theta = mPart->Theta();
 825         if  (fDebug>=2 && mPart->GetStatusCode()==1) {
 826            printf("ind %7i part %7i -> child1 %5i child2 %5i Status %5i\n",
 827            part, mpart, child1, mPart->GetLastDaughter(), mPart->GetStatusCode());
 828         }
 829
 830         if (fDebug >= 2 && genType == kPythia) {
 831             if (part == 6 || part == 7)
 832             {
 833                 printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f\n",
 834                        part-5, pT, eta, phi);
 835             }
 836         }
 837
 838         fTrackList[part] = 0;
 839         fPtT[part]       = pT; // must be change after correction for resolution !!!
 840         fEtaT[part]      = eta;
 841         fPhiT[part]      = phi;
 842         fPdgT[part]      = mpart;
 843
 844 //      final state particles only
 845
 846         if (genType == kPythia) {
 847             if (mPart->GetStatusCode() != 1) continue;
 848         } else if (genType == kHijing) {
 849             if (child1 >= 0 && child1 < npart) continue;
 850         }
 851
 852
 853         TParticlePDG* pdgP = 0;
 854 // charged or neutral
 855         pdgP  = mPart->GetPDG();
 856         chTmp = pdgP->Charge() / 3.; // 13-feb-2001!!
 857
 858         if (ich == 0) {
 859             if (chTmp == 0) {
 860                 if (!fK0N) {
 861                     continue;
 862                 } else {
 863                     if (mpart != kNeutron    &&
 864                         mpart != kNeutronBar &&
 865                         mpart != kK0Long) continue;
 866                 }
 867             }
 868         } else if (ich == 2) {
 869           if (mpart == kNeutron    ||
 870               mpart == kNeutronBar ||
 871               mpart == kK0Long) continue;
 872         }
 873
 874         if (TMath::Abs(eta)<=0.9) fNChTpc++;
 875 // final state only
 876         if (child1 >= 0 && child1 < npart) continue;
 877 // acceptance cut
 878         if (eta > fEtaMax || eta < fEtaMin)    continue;
 879         if (phi > fPhiMax || phi < fPhiMin )   continue;
 880 //
 881         if (fDebug >= 3)
 882         printf("\n=>nsel:%5d mpart %5d child1 %5d eta %6.2f phi %6.2f pT %6.2f ch %3.0f ",
 883         part, mpart, child1, eta, phi, pT, chTmp);
 884 //
 885 //
 886 // Momentum smearing goes here ...
 887 //
 888         fhTrackPt->Fill(pT);
 889         Float_t pw;
 890         if (fSmear && TMath::Abs(chTmp)) {
 891             pw = AliEMCALFast::SmearMomentum(1,p);
 892             // p changed - take into account when calculate pT,
 893             // pz and so on ..  ?
 894             pT = (pw/p) * pT;
 895             if(fDebug >= 4) printf("\n Smearing : p %8.4f change to %8.4f ", p, pw);
 896             p  = pw;
 897         }
 898 //
 899 // Tracking Efficiency and TPC acceptance goes here ...
 900         Float_t eff;
 901         if (fEffic && TMath::Abs(chTmp)) {
 902             eff = 0.9; //  AliEMCALFast::Efficiency(2,p);
 903             if(fhEff) fhEff->Fill(p, eff);
 904             if (AliEMCALFast::RandomReject(eff)) {
 905                 if(fDebug >= 5) printf(" reject due to unefficiency ");
 906                 continue;
 907             }
 908         }
 909 //
 910 // Correction of Hadronic Energy goes here ...
 911 //
 912 //
 913 // phi propagation for hadronic correction
 914
 915         Bool_t curls = kFALSE; // hit two the EMCAL (no curl)
 916         Float_t phiHC=0.0, dpH=0.0, dphi=0.0, eTdpH=0;
 917         if(TMath::Abs(chTmp)) {
 918             // hadr. correction only for charge particle
 919             dphi  = PropagatePhi(pT, chTmp, curls);
 920             phiHC = phi + dphi;
 921             if (fDebug >= 6) {
 922                 printf("\n Delta phi %f pT %f ", dphi, pT);
 923                 if (curls) printf("\n !! Track is curling");
 924             }
 925             if(!curls) fhTrackPtBcut->Fill(pT);
 926
 927             if (fHCorrection && !curls) {
 928                 if (!fHadronCorrector)
 929                     Fatal("AliEMCALJetFinder",
 930                           "Hadronic energy correction required but not defined !");
 931
 932                 dpH    = fHadronCorrector->GetEnergy(p, eta, 7);
 933                 eTdpH  = dpH*TMath::Sin(theta);
 934
 935                 if (fDebug >= 7) printf(" phi %f phiHC %f eTcorr %f\n",
 936                                         phi, phiHC, -eTdpH); // correction is negative
 937                 Int_t xbin,ybin;
 938                 xbin = fLego->GetXaxis()->FindBin(eta);
 939                 ybin = fLego->GetYaxis()->FindBin(phiHC);
 940                 cout <<"Hadron Correction affected bin - contents before correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
 941                 fLego->Fill(eta, phi, -fSamplingF*eTdpH );
 942                 cout <<"Hadron Correction affected bin - contents after  correction : "<<fLego->GetBinContent(xbin,ybin)<<endl;
 943                 fhLegoHadrCorr->Fill(eta, phi, fSamplingF*eTdpH);
 944             }
 945         }
 946 //
 947 //  More to do ? Just think about it !
 948 //
 949         if (phi > fPhiMax || phi < fPhiMin )   continue;
 950
 951         if(TMath::Abs(chTmp) ) { // charge particle
 952             if (pT > fPtCut && !curls) {
 953                 if (fDebug >= 8) printf("Charge :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
 954                                         eta , phi, pT, fNtS);
 955                 fLego->Fill(eta, phi, pT);
 956                 fhLegoTracks->Fill(eta, phi, pT); // 20-feb for checking
 957                 fTrackList[part] = 1;
 958                 fNtS++;
 959             }
 960         } else if(ich > 0 || fK0N) {
 961             // case of n, nbar and K0L
 962             if (fDebug >= 9) printf("Neutral :  fLego->Fill(%5.2f, %5.2f, %6.2f, %d)\n",
 963                                     eta , phi, pT, fNtS);
 964             fLego->Fill(eta, phi, pT);
 965             fTrackList[part] = 1;
 966             fNtS++;
 967         }
 968     } // primary loop
 969     Float_t etsum = 0.;
 970     for(Int_t i=0; i<fLego->GetSize(); i++) {
 971         Float_t etc =  (*fLego)[i];
 972         if (etc > fMinCellEt) etsum += etc;
 973     }
 974
 975     printf("FillFromTracks: Sum above threshold %f -> %f (%f)\n", fMinCellEt, etsum, fLego->Integral());
 976     printf("                Track selected(fNtS) %i \n", fNtS);
 977
 978     DumpLego();
 979 }
 980
 981 void AliEMCALJetFinder::FillFromHits(Int_t flag)
 982 {
 983 //
 984 // Fill Cells with hit information
 985 //
 986 //
 987     if (fDebug >= 2)
 988     printf("\n AliEMCALJetFinder::FillFromHits(%i)\n",flag);
 989
 990     ResetMap();
 991
 992     if (!fLego) BookLego();
 993 //  Reset eta-phi maps if needed
 994     if (flag == 0) { // default behavior
 995       fLego->Reset();
 996       fhLegoTracks->Reset();
 997       fhLegoEMCAL->Reset();
 998       fhLegoHadrCorr->Reset();
 999     }
1000 //  Initialize from background event if available
1001 //
1002 // Access hit information
1003     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1004     AliLoader *emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL");
1005     TTree *treeH = emcalLoader->TreeH();
1006     Int_t ntracks = (Int_t) treeH->GetEntries();
1007 //
1008 //   Loop over tracks
1009 //
1010     Int_t nbytes = 0;
1011     //    Double_t etH = 0.0;
1012
1013     for (Int_t track=0; track<ntracks;track++) {
1014         gAlice->ResetHits();
1015         nbytes += treeH->GetEvent(track);
1016 //
1017 //  Loop over hits
1018 //
1019         for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1020             mHit;
1021             mHit=(AliEMCALHit*) pEMCAL->NextHit())
1022         {
1023             Float_t x      =    mHit->X();         // x-pos of hit
1024             Float_t y      =    mHit->Y();         // y-pos
1025             Float_t z      =    mHit->Z();         // z-pos
1026             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
1027
1028             Float_t r      =    TMath::Sqrt(x*x+y*y);
1029             Float_t theta  =    TMath::ATan2(r,z);
1030             Float_t eta    =   -TMath::Log(TMath::Tan(theta/2.));
1031             Float_t phi    =    TMath::ATan2(y,x);
1032
1033             if (fDebug >= 21) printf("\n Hit %f %f %f %f %f %f %f %f", x, y, z, eloss, r, eta, phi, fSamplingF);
1034
1035             //            etH = fSamplingF*eloss*TMath::Sin(theta);
1036             fLego->Fill(eta, phi, eloss);
1037         } // Hit Loop
1038     } // Track Loop
1039
1040     // Transition from deposit energy to eT (eT = de*SF*sin(theta))
1041     Double_t etsum = 0;
1042     for(Int_t i=1; i<=fLego->GetNbinsX(); i++){   // eta
1043        Double_t sinTheta = fhSinTheta->GetBinContent(i), eT=0;
1044        for(Int_t j=1; j<=fLego->GetNbinsY(); j++){ // phi
1045           eT = fLego->GetBinContent(i,j)*fSamplingF*sinTheta;
1046           fLego->SetBinContent(i,j,eT);
1047     // copy content of fLego to fhLegoEMCAL (fLego and fhLegoEMCAL are identical)
1048           fhLegoEMCAL->SetBinContent(i,j,eT);
1049           if (eT > fMinCellEt) etsum += eT;
1050        }
1051     }
1052
1053     //    for(Int_t i=0; i<fLego->GetSize(); i++) {
1054     //  (*fhLegoEMCAL)[i] = (*fLego)[i];
1055     //  Float_t etc =  (*fLego)[i];
1056     //  if (etc > fMinCellEt) etsum += etc;
1057     //    }
1058
1059     printf("FillFromHits: Sum above threshold %f -> %f \n ", fMinCellEt, etsum);
1060     //    DumpLego(); ??
1061 }
1062
1063 void AliEMCALJetFinder::FillFromDigits(Int_t flag)
1064 {
1065 //
1066 // Fill Cells with digit information
1067 //
1068 //
1069     ResetMap();
1070
1071     if (!fLego) BookLego();
1072     if (flag == 0) fLego->Reset();
1073     Int_t nbytes=0;
1074
1075
1076 //
1077 //  Connect digits
1078 //
1079     TClonesArray* digs      = new TClonesArray("AliEMCALDigit", 10000);
1080     TTree *treeD = gAlice->TreeD();
1081     TBranchElement* branchDg = (TBranchElement*)
1082         treeD->GetBranch("EMCAL");
1083
1084     if (!branchDg) Fatal("AliEMCALJetFinder",
1085                          "Reading digits requested but no digits in file !");
1086
1087     branchDg->SetAddress(&digs);
1088     Int_t nent = (Int_t) branchDg->GetEntries();
1089 //
1090 //  Connect digitizer
1091 //
1092     AliEMCALDigitizer* digr = new AliEMCALDigitizer();
1093     TBranchElement* branchDr = (TBranchElement*)
1094         treeD->GetBranch("AliEMCALDigitizer");
1095     branchDr->SetAddress(&digr);
1096 //
1097 //
1098     nbytes += branchDg->GetEntry(0);
1099     nbytes += branchDr->GetEntry(0);
1100 //
1101 //  Get digitizer parameters
1102     Float_t ecADCped  = digr->GetECApedestal();
1103     Float_t ecADCcha  = digr->GetECAchannel();
1104
1105     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1106     AliEMCALGeometry* geom =
1107         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
1108
1109     if (fDebug) {
1110         Int_t   ndig = digs->GetEntries();
1111         Info("FillFromDigits","Number of Digits: %d %d\n Parameters: EC: %f %f\n Geometry: %d %d",
1112              ndig, nent, ecADCped, ecADCcha, geom->GetNEta(), geom->GetNPhi());
1113     }
1114
1115 //
1116 //  Loop over digits
1117     AliEMCALDigit* sdg;
1118     TIter next(digs);
1119     while ((sdg = (AliEMCALDigit*)(next())))
1120     {
1121         Double_t pedestal = 0.;
1122         Double_t channel  = 0.;
1123         if (geom->CheckAbsCellId(sdg->GetId()))  { // May 31, 2006; IsInECA(Int_t) -> CheckAbsCellId
1124           pedestal = ecADCped;
1125           channel  = ecADCcha;
1126         }
1127         else
1128           Fatal("FillFromDigits", "unexpected digit is number!") ;
1129
1130         Float_t eta = sdg->GetEta();
1131         Float_t phi = sdg->GetPhi() * TMath::Pi()/180.;
1132         Float_t amp = (Float_t) (channel*(sdg->GetAmp())-pedestal);
1133
1134         if (fDebug > 1)
1135           Info("FillFromDigits", "Digit: eta %8.3f phi %8.3f amp %8.3f %8d",
1136                            eta, phi, amp, sdg->GetAmp());
1137
1138         fLego->Fill(eta, phi, fSamplingF*amp);
1139     } // digit loop
1140 //
1141 //  Dump lego hist
1142     DumpLego();
1143 }
1144
1145
1146 void AliEMCALJetFinder::FillFromHitFlaggedTracks(Int_t flag)
1147 {
1148 //
1149 // Fill Cells with hit information
1150 //
1151 //
1152     ResetMap();
1153
1154     if (!fLego) BookLego();
1155 // Reset
1156     if (flag == 0) fLego->Reset();
1157
1158 // Flag charged tracks passing through TPC which
1159 // are associated to EMCAL Hits
1160     BuildTrackFlagTable();
1161
1162 //
1163 // Access particle information
1164     TTree *treeK = gAlice->TreeK();
1165     Int_t ntracks = (Int_t) treeK->GetEntries();
1166
1167     if (fPtT)       delete[] fPtT;
1168     if (fEtaT)      delete[] fEtaT;
1169     if (fPhiT)      delete[] fPhiT;
1170     if (fPdgT)      delete[] fPdgT;
1171
1172     fPtT       = new Float_t[ntracks];
1173     fEtaT      = new Float_t[ntracks];
1174     fPhiT      = new Float_t[ntracks];
1175     fPdgT      = new Int_t[ntracks];
1176
1177     fNt   = ntracks;
1178     fNtS  = 0;
1179
1180     for (Int_t track = 0; track < ntracks; track++) {
1181           TParticle *mPart = gAlice->GetMCApp()->Particle(track);
1182           Float_t pT    = mPart->Pt();
1183           Float_t phi   = mPart->Phi();
1184           Float_t eta   = mPart->Eta();
1185
1186         if(fTrackList[track]) {
1187           fPtT[track]       = pT;
1188           fEtaT[track]      = eta;
1189           fPhiT[track]      = phi;
1190           fPdgT[track]      = mPart->GetPdgCode();
1191
1192           if (track < 2) continue;      //Colliding particles?
1193           if (pT == 0 || pT < fPtCut) continue;
1194           fNtS++;
1195           fLego->Fill(eta, phi, pT);
1196         }
1197       } // track loop
1198     DumpLego();
1199 }
1200
1201 void AliEMCALJetFinder::FillFromParticles()
1202 {
1203 // 26-feb-2002 PAI - for checking all chain
1204 // Work on particles level; accept all particle (not neutrino )
1205
1206     Double_t pX=0, pY=0, pZ=0, energy=0; // checking conservation law
1207     fNChTpc = 0;
1208
1209     ResetMap();
1210     if (!fLego) BookLego();
1211     fLego->Reset();
1212 //
1213 // Access particles information
1214     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1215     if (fDebug >= 2 || npart<=0) {
1216        printf(" AliEMCALJetFinder::FillFromParticles : npart %i\n", npart);
1217        if(npart<=0) return;
1218     }
1219     fNt   = npart;
1220     fNtS  = 0;
1221     RearrangeParticlesMemory(npart);
1222
1223 //  Go through the particles
1224     Int_t mpart, child1, child2, geantPdg;
1225     Float_t pT, phi, eta, e=0, px=0, py=0, pz=0;
1226     TParticle *mPart=0;
1227     for (Int_t part = 0; part < npart; part++) {
1228
1229         fTrackList[part] = 0;
1230
1231         mPart   = gAlice->GetMCApp()->Particle(part);
1232         mpart   = mPart->GetPdgCode();
1233         child1  = mPart->GetFirstDaughter();
1234         child2  = mPart->GetLastDaughter();
1235         pT      = mPart->Pt();
1236         phi     = mPart->Phi();
1237         eta     = mPart->Eta();
1238
1239         px      = mPart->Px();
1240         py      = mPart->Py();
1241         pz      = mPart->Pz();
1242         e       = mPart->Energy();
1243
1244 // see pyedit in Pythia's text
1245         geantPdg = mpart;
1246         if (IsThisPartonsOrDiQuark(mpart)) continue;
1247         printf("%5i: %5i(%2i) px %5.1f py %5.1f pz %6.1f e %6.1f childs %5i,%5i \n",
1248         part, mpart, geantPdg, px, py, pz, e, child1, child2);
1249
1250 //  exclude partons (21 - gluon, 92 - string)
1251
1252
1253 // exclude neutrinous also ??
1254         if (fDebug >= 11 && pT>0.01)
1255         printf("\n part:%5d mpart %5d eta %9.2f phi %9.2f pT %9.2f ",
1256         part, mpart, eta, phi, pT);
1257
1258         fPtT[part]       = pT;
1259         fEtaT[part]      = eta;
1260         fPhiT[part]      = phi;
1261         fPdgT[part]      = mpart;
1262         fNtS++;
1263
1264 // final state only
1265         if (child1 >= 0 && child1 < npart) continue;
1266
1267         //        printf("%4i -> %5i(%3i) px %6.1f py %6.1f pz %7.1f e %8.2f child1 %5i %s\n",
1268         //        part, mpart, geantPdg, px, py, pz, e, child1, name.Data());
1269         pX += px;
1270         pY += py;
1271         pZ += pz;
1272         energy  += e;
1273
1274
1275         if (TMath::Abs(eta) <= 0.9) fNChTpc++;
1276 // acceptance cut
1277         if (eta > fEtaMax || eta < fEtaMin)    continue;
1278         if (phi > fPhiMax || phi < fPhiMin )   continue;
1279 //
1280         if(fK0N==0 ) { // exclude neutral hadrons
1281           if (mpart == kNeutron || mpart == kNeutronBar || mpart == kK0Long) continue;
1282         }
1283         fTrackList[part] = 1;
1284         fLego->Fill(eta, phi, pT);
1285
1286     } // primary loop
1287     printf("\n                PX %8.2f  PY %8.2f  PZ %8.2f  E %8.2f \n",
1288     pX, pY, pZ, energy);
1289     DumpLego();
1290     if(fhChPartMultInTpc) fhChPartMultInTpc->Fill(fNChTpc);
1291 }
1292
1293 void AliEMCALJetFinder::FillFromPartons()
1294 {
1295 // function under construction - 13-feb-2002 PAI
1296
1297     if (fDebug >= 2)
1298     printf("\n AliEMCALJetFinder::FillFromPartons()\n");
1299     //
1300
1301     ResetMap();
1302     if (!fLego) BookLego();
1303     fLego->Reset();
1304 //
1305 // Access particle information
1306     Int_t npart = (gAlice->GetHeader())->GetNprimary();
1307     if (fDebug >= 2 || npart<=0)
1308     printf("\n AliEMCALJetFinder::FillFromPartons : npart %i\n", npart);
1309     fNt   = 0; // for FindTracksInJetCone
1310     fNtS  = 0;
1311
1312 //  Go through the partons
1313     Int_t statusCode=0;
1314     for (Int_t part = 8; part < npart; part++) {
1315         TParticle *mPart = gAlice->GetMCApp()->Particle(part);
1316         Int_t mpart   = mPart->GetPdgCode();
1317         //      Int_t child1  = MPart->GetFirstDaughter();
1318         Float_t pT    = mPart->Pt();
1319         //      Float_t p     = MPart->P();
1320         Float_t phi   = mPart->Phi();
1321         Float_t eta   = mPart->Eta();
1322         //      Float_t theta = MPart->Theta();
1323         statusCode    = mPart->GetStatusCode();
1324
1325 // accept partons (21 - gluon, 92 - string)
1326         if (!(TMath::Abs(mpart) <= 6 || mpart == 21 ||mpart == 92)) continue;
1327         if (fDebug > 1 && pT>0.01)
1328         printf("\n part:%5d mpart %5d status  %5d eta %8.2f phi %8.2f pT %8.2f ",
1329         part, mpart, statusCode, eta, phi, pT);
1330         //      if (fDebug >= 3) MPart->Print();
1331 // accept partons before fragmentation - p.57 in Pythia manual
1332 //        if(statusCode != 1) continue;
1333 // acceptance cut
1334         if (eta > fEtaMax || eta < fEtaMin)    continue;
1335         if (phi > fPhiMax || phi < fPhiMin )   continue;
1336 // final state only
1337 //      if (child1 >= 0 && child1 < npart) continue;
1338 //
1339 //
1340         fLego->Fill(eta, phi, pT);
1341
1342     } // primary loop
1343     DumpLego();
1344 }
1345
1346 void AliEMCALJetFinder::BuildTrackFlagTable() {
1347
1348 // Method to generate a lookup table for TreeK particles
1349 // which are linked to hits in the EMCAL
1350 //
1351 // --Author: J.L. Klay
1352 //
1353 // Access hit information
1354     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
1355
1356     TTree *treeK = gAlice->TreeK();             // Get the number of entries in the kine tree
1357     Int_t nKTrks = (Int_t) treeK->GetEntries(); // (Number of particles created somewhere)
1358
1359     if(fTrackList) delete[] fTrackList;         //Make sure we get rid of the old one
1360     fTrackList = new Int_t[nKTrks];             //before generating a new one
1361
1362     for (Int_t i = 0; i < nKTrks; i++) {        //Initialize members to 0
1363         fTrackList[i] = 0;
1364     }
1365
1366     AliLoader *emcalLoader = AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL");
1367     //AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
1368     // TTree *treeH = gAlice->TreeH();
1369     TTree *treeH = emcalLoader->TreeH();
1370    Int_t ntracks = (Int_t) treeH->GetEntries();
1371 //
1372 //   Loop over tracks
1373 //
1374     Int_t nbytes = 0;
1375
1376     for (Int_t track=0; track<ntracks;track++) {
1377         gAlice->ResetHits();
1378         nbytes += treeH->GetEvent(track);
1379         if (pEMCAL)  {
1380
1381 //
1382 //  Loop over hits
1383 //
1384             for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
1385                 mHit;
1386                 mHit=(AliEMCALHit*) pEMCAL->NextHit())
1387             {
1388                 Int_t   iTrk   =    mHit->Track();       // track number
1389                 Int_t   idprim =    mHit->GetPrimary();  // primary particle
1390
1391                 //Determine the origin point of this particle - it made a hit in the EMCAL
1392                 TParticle         *trkPart = gAlice->GetMCApp()->Particle(iTrk);
1393                 TParticlePDG  *trkPdg  = trkPart->GetPDG();
1394                 Int_t      trkCode = trkPart->GetPdgCode();
1395                 Double_t           trkChg;
1396                 if (trkCode < 10000) {          //Big Ions cause problems for
1397                     trkChg  = trkPdg->Charge(); //this function.  Since they aren't
1398                 } else {                                //likely to make it very far, set
1399                     trkChg  = 0.0;                      //their charge to 0 for the Flag test
1400                 }
1401                 Float_t            vxTrk   = trkPart->Vx();
1402                 Float_t            vyTrk   = trkPart->Vy();
1403                 Float_t            vrTrk   = TMath::Sqrt(vxTrk*vxTrk+vyTrk*vyTrk);
1404                 fTrackList[iTrk] = SetTrackFlag(vrTrk,trkCode,trkChg);
1405
1406                 //Loop through the ancestry of the EMCAL entrance particles
1407                 Int_t ancestor = trkPart->GetFirstMother();  //Get track's Mother
1408                 while (ancestor != -1) {
1409                     TParticle     *ancPart = gAlice->GetMCApp()->Particle(ancestor);  //get particle info on ancestor
1410                     TParticlePDG  *ancPdg  = ancPart->GetPDG();
1411                     Int_t             ancCode = ancPart->GetPdgCode();
1412                     Double_t       ancChg;
1413                     if (ancCode < 10000) {
1414                         ancChg  = ancPdg->Charge();
1415                     } else {
1416                         ancChg  = 0.0;
1417                     }
1418                     Float_t           vxAnc   = ancPart->Vx();
1419                     Float_t           vyAnc   = ancPart->Vy();
1420                     Float_t           vrAnc   = TMath::Sqrt(vxAnc*vxAnc+vyAnc*vyAnc);
1421                     fTrackList[ancestor] = SetTrackFlag(vrAnc,ancCode,ancChg);
1422                     ancestor = ancPart->GetFirstMother();           //Get the next ancestor
1423                 }
1424
1425                 //Determine the origin point of the primary particle associated with the hit
1426                 TParticle     *primPart = gAlice->GetMCApp()->Particle(idprim);
1427                 TParticlePDG  *primPdg  = primPart->GetPDG();
1428                 Int_t      primCode = primPart->GetPdgCode();
1429                 Double_t           primChg;
1430                 if (primCode < 10000) {
1431                     primChg  = primPdg->Charge();
1432                 } else {
1433                     primChg  = 0.0;
1434                 }
1435                 Float_t    vxPrim = primPart->Vx();
1436                 Float_t    vyPrim = primPart->Vy();
1437                 Float_t    vrPrim = TMath::Sqrt(vxPrim*vxPrim+vyPrim*vyPrim);
1438                 fTrackList[idprim] = SetTrackFlag(vrPrim,primCode,primChg);
1439             } // Hit Loop
1440         } //if (pEMCAL)
1441     } // Track Loop
1442 }
1443
1444 Int_t AliEMCALJetFinder
1445 ::SetTrackFlag(Float_t radius, Int_t code, Double_t charge) {
1446 // Set the flag for the track
1447     Int_t flag    = 0;
1448     Int_t parton  = 0;
1449     Int_t neutral = 0;
1450
1451     if (charge == 0) neutral = 1;
1452
1453     if (TMath::Abs(code) <= 6   ||
1454         code == 21              ||
1455         code == 92) parton = 1;
1456
1457     //It's not a parton, it's charged and it went through the TPC
1458     if (!parton && !neutral && radius <= 84.0) flag = 1;
1459
1460     return flag;
1461 }
1462
1463
1464
1465 void AliEMCALJetFinder::SaveBackgroundEvent(const char *name)
1466 {
1467 // Saves the eta-phi lego and the tracklist and name of file with BG events
1468 //
1469     if (fLegoB) {
1470        fLegoB->Reset();
1471        (*fLegoB) = (*fLegoB) + (*fLego);
1472        //       if(fDebug)
1473        printf("\n AliEMCALJetFinder::SaveBackgroundEvent() (fLegoB) %f = %f(fLego) \n",
1474        fLegoB->Integral(), fLego->Integral());
1475     }
1476
1477     if (fPtB)        delete[] fPtB;
1478     if (fEtaB)       delete[] fEtaB;
1479     if (fPhiB)       delete[] fPhiB;
1480     if (fPdgB)       delete[] fPdgB;
1481     if (fTrackListB) delete[] fTrackListB;
1482
1483     fPtB          = new Float_t[fNtS];
1484     fEtaB         = new Float_t[fNtS];
1485     fPhiB         = new Float_t[fNtS];
1486     fPdgB         = new Int_t  [fNtS];
1487     fTrackListB   = new Int_t  [fNtS];
1488
1489     fNtB = 0;
1490
1491     for (Int_t i = 0; i < fNt; i++) {
1492         if (!fTrackList[i]) continue;
1493         fPtB [fNtB]       = fPtT [i];
1494         fEtaB[fNtB]       = fEtaT[i];
1495         fPhiB[fNtB]       = fPhiT[i];
1496         fPdgB[fNtB]       = fPdgT[i];
1497         fTrackListB[fNtB] = 1;
1498         fNtB++;
1499     }
1500     fBackground = 1;
1501     printf(" fNtB %i => fNtS %i #particles %i \n", fNtB, fNtS, fNt);
1502
1503     if(strlen(name) == 0) {
1504        TSeqCollection *li = gROOT->GetListOfFiles();
1505        TString nf;
1506        for(Int_t i=0; i<li->GetSize(); i++) {
1507           nf = ((TFile*)li->At(i))->GetName();
1508           if(nf.Contains("backgorund")) break;
1509        }
1510        fBGFileName = nf;
1511     } else {
1512        fBGFileName = name;
1513     }
1514     printf("BG file name is \n %s\n", fBGFileName.Data());
1515 }
1516
1517 void AliEMCALJetFinder::InitFromBackground()
1518 {
1519 //
1520 //
1521     if (fDebug) printf("\n AliEMCALJetFinder::InitFromBackground() ");
1522
1523     if (fLego) {
1524         fLego->Reset();
1525         (*fLego) = (*fLego) + (*fLegoB);
1526         if(fDebug)
1527             printf("\n AliEMCALJetFinder::InitBackgroundEvent() (fLego) %f = %f(fLegoB) \n",
1528                    fLego->Integral(), fLegoB->Integral());
1529     } else {
1530         printf(" => fLego undefined \n");
1531     }
1532 }
1533
1534
1535 void AliEMCALJetFinder::FindTracksInJetCone()
1536 {
1537 //
1538 //  Build list of tracks inside jet cone
1539 //
1540 //  Loop over jets
1541     Int_t njet = Njets();
1542     for (Int_t nj = 0; nj < njet; nj++)
1543     {
1544         Float_t etaj = JetEtaW(nj);
1545         Float_t phij = JetPhiW(nj);
1546         Int_t   nT   = 0;           // number of associated tracks
1547
1548 // Loop over particles in current event
1549         for (Int_t part = 0; part < fNt; part++) {
1550             if (!fTrackList[part]) continue;
1551             Float_t phi      = fPhiT[part];
1552             Float_t eta      = fEtaT[part];
1553             Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1554                                            (phij-phi)*(phij-phi));
1555             if (dr < fConeRadius) {
1556                 fTrackList[part] = nj+2;
1557                 nT++;
1558             } // < ConeRadius ?
1559         } // particle loop
1560
1561 // Same for background event if available
1562         Int_t nTB = 0;
1563         if (fBackground) {
1564             for (Int_t part = 0; part < fNtB; part++) {
1565                 Float_t phi      = fPhiB[part];
1566                 Float_t eta      = fEtaB[part];
1567                 Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
1568                                                (phij-phi)*(phij-phi));
1569                 fTrackListB[part] = 1;
1570
1571                 if (dr < fConeRadius) {
1572                     fTrackListB[part] = nj+2;
1573                     nTB++;
1574                 } // < ConeRadius ?
1575             } // particle loop
1576         } // Background available ?
1577
1578         Int_t nT0 = nT + nTB;
1579         printf("Total number of tracks %d\n", nT0);
1580
1581         if (nT0 > 1000) nT0 = 1000;
1582
1583         Float_t* ptT  = new Float_t[nT0];
1584         Float_t* etaT = new Float_t[nT0];
1585         Float_t* phiT = new Float_t[nT0];
1586         Int_t*   pdgT = new Int_t[nT0];
1587
1588         Int_t iT = 0;
1589         Int_t j;
1590
1591         for (Int_t part = 0; part < fNt; part++) {
1592             if (fTrackList[part] == nj+2) {
1593                 Int_t index = 0;
1594                 for (j=iT-1; j>=0; j--) {
1595                     if (fPtT[part] > ptT[j]) {
1596                         index = j+1;
1597                         break;
1598                     }
1599                 }
1600                 for (j=iT-1; j>=index; j--) {
1601                     ptT [j+1]  = ptT [j];
1602                     etaT[j+1]  = etaT[j];
1603                     phiT[j+1]  = phiT[j];
1604                     pdgT[j+1]  = pdgT[j];
1605                 }
1606                 ptT [index] = fPtT [part];
1607                 etaT[index] = fEtaT[part];
1608                 phiT[index] = fPhiT[part];
1609                 pdgT[index] = fPdgT[part];
1610                 iT++;
1611             } // particle associated
1612             if (iT > nT0) break;
1613         } // particle loop
1614
1615         if (fBackground) {
1616             for (Int_t part = 0; part < fNtB; part++) {
1617                 if (fTrackListB[part] == nj+2) {
1618                     Int_t index = 0;
1619                     for (j=iT-1; j>=0; j--) {
1620                         if (fPtB[part] > ptT[j]) {
1621                             index = j+1;
1622
1623                             break;
1624                         }
1625                     }
1626                     for (j=iT-1; j>=index; j--) {
1627                         ptT [j+1]  = ptT [j];
1628                         etaT[j+1]  = etaT[j];
1629                         phiT[j+1]  = phiT[j];
1630                         pdgT[j+1]  = pdgT[j];
1631                     }
1632                     ptT [index] = fPtB [part];
1633                     etaT[index] = fEtaB[part];
1634                     phiT[index] = fPhiB[part];
1635                     pdgT[index] = fPdgB[part];
1636                     iT++;
1637                 } // particle associated
1638                 if (iT > nT0) break;
1639             } // particle loop
1640         } // Background available ?
1641
1642         fJetT[nj]->SetTrackList(nT0, ptT, etaT, phiT, pdgT);
1643         delete[] ptT;
1644         delete[] etaT;
1645         delete[] phiT;
1646         delete[] pdgT;
1647
1648     } // jet loop loop
1649 }
1650
1651 Float_t AliEMCALJetFinder::PropagatePhi(Float_t pt, Float_t charge, Bool_t& curls)
1652 {
1653 // Propagates phi angle to EMCAL radius
1654 //
1655   static Float_t b = 0.0, rEMCAL = -1.0;
1656 // Get field in kGS
1657   b = gAlice->Field()->SolenoidField();
1658 // Get EMCAL radius in cm
1659   rEMCAL = AliEMCALGeometry::GetInstance()->GetIPDistance();
1660   Float_t dPhi = 0.;
1661 //
1662 //
1663 // bending radies
1664 // pt [Gev]
1665 // B  [kG]
1666 //
1667   Float_t rB = 3335.6 * pt / TMath::Abs(b);  // [cm]  (case of |charge|=1)
1668 //
1669 // check if particle is curling below EMCAL
1670     if (2.*rB < rEMCAL) {
1671         curls = kTRUE;
1672         AliDebug(1, Form(" Low pt %f \n", pt));
1673         return dPhi;
1674     }
1675 //
1676 // if not calculate delta phi
1677     Float_t phi = TMath::ACos(1.-rEMCAL*rEMCAL/(2.*rB*rB));
1678     dPhi = TMath::ATan2(1.-TMath::Cos(phi), TMath::Sin(phi));
1679     dPhi = -TMath::Sign(dPhi, charge);
1680 //
1681     AliDebug(1, Form(" B %7.3f kG : rEMCAL %7.2f : dphi %7.5f(%7.5f)\n",
1682                      b, rEMCAL, dPhi,dPhi*TMath::RadToDeg()));
1683     return dPhi;
1684 }
1685
1686 void hf1(Int_t& , Float_t& , Float_t& )
1687 {
1688 // dummy for hbook calls
1689     ;
1690 }
1691
1692 void AliEMCALJetFinder::DrawLego(const char *opt)
1693 {
1694 // Draw lego
1695         if(fLego) fLego->Draw(opt);
1696 }
1697
1698 void  AliEMCALJetFinder::DrawLegoBackground(const char *opt)
1699 {
1700 // Draw background lego
1701         if(fLegoB) fLegoB->Draw(opt);
1702 }
1703
1704 void AliEMCALJetFinder::DrawLegoEMCAL(const char *opt)
1705 {
1706 // Draw EMCAL Lego
1707         if(fhLegoEMCAL) fhLegoEMCAL->Draw(opt);
1708 }
1709
1710 void AliEMCALJetFinder::DrawHistsForTuning(Int_t mode)
1711 {
1712 // Draw all hists
1713   static TPaveText *varLabel=0;
1714   if(!fC1) {
1715     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1716   }
1717   fC1->Clear();
1718   fC1->Divide(2,2);
1719   fC1->cd(1);
1720   gPad->SetLogy(1);
1721   fhCellEt->Draw();
1722
1723   fC1->cd(2);
1724   gPad->SetLogy(1);
1725   fhCellEMCALEt->Draw();
1726
1727   fC1->cd(3);
1728   gPad->SetLogy(1);
1729   fhTrackPt->Draw();
1730   fhTrackPtBcut->SetLineColor(2);
1731   fhTrackPtBcut->Draw("same");
1732
1733   fC1->cd(4);
1734   if(!varLabel) {
1735     PrintParameters(1);
1736     varLabel = new TPaveText(0.05,0.5,0.95,0.95,"NDC");
1737     varLabel->SetTextAlign(12);
1738     varLabel->SetFillColor(19); // see TAttFill
1739     TString tmp(GetTitle());
1740     varLabel->ReadFile(GetFileNameForParameters());
1741   }
1742   varLabel->Draw();
1743   fC1->Update();
1744   if(mode) { // for saving picture to the file
1745     TString stmp(GetFileNameForParameters());
1746     stmp.ReplaceAll("_Par.txt",".ps");
1747     fC1->Print(stmp.Data());
1748   }
1749 }
1750
1751 void AliEMCALJetFinder::PrintParameters(Int_t mode)
1752 {
1753 // Print all parameters out
1754
1755   FILE *file=0;
1756   if(mode==0) file = stdout; // output to terminal
1757   else {
1758     file = fopen(GetFileNameForParameters(),"w");
1759     if(file==0) file = stdout;
1760   }
1761   fprintf(file,"====   Filling lego   ==== \n");
1762   fprintf(file,"Sampling fraction %6.3f  ", fSamplingF);
1763   fprintf(file,"Smearing          %6i  ", fSmear);
1764   fprintf(file,"Efficiency        %6i\n", fEffic);
1765   fprintf(file,"Hadr.Correct.     %6i  ", fHCorrection);
1766   fprintf(file,"P_{T} Cut of ch.par. %6.2f\n", fPtCut);
1767   fprintf(file,"====  Jet finding     ==== \n");
1768   fprintf(file,"Cone radius       %6.2f  ", fConeRadius);
1769   fprintf(file,"Seed E_{T}           %6.1f\n", fEtSeed);
1770   fprintf(file,"Min E_{T} of cell    %6.1f  ", fMinCellEt);
1771   fprintf(file,"Min E_{T} of jet     %6.1f\n", fMinJetEt);
1772   if(fMode) {
1773     fprintf(file,"====  Bg subtraction     ==== \n");
1774     fprintf(file,"BG subtraction  %6i  ", fMode);
1775     fprintf(file,"Min cone move   %6.2f\n", fMinMove);
1776     fprintf(file,"Max cone move   %6.2f  ", fMaxMove);
1777     fprintf(file,"%% change for BG %6.4f\n", fPrecBg);
1778   } else
1779   fprintf(file,"==== No Bg subtraction     ==== \n");
1780   //
1781   printf("BG file name is %s \n", fBGFileName.Data());
1782   if(file != stdout) fclose(file);
1783 }
1784
1785 void AliEMCALJetFinder::DrawLegos()
1786 {
1787 // Draw all legos
1788   if(!fC1) {
1789     fC1 = new TCanvas("C1","Hists. for tunning", 0,25,600,800);
1790   }
1791   fC1->Clear();
1792   fC1->Divide(2,2);
1793   gStyle->SetOptStat(111111);
1794
1795   Int_t nent1, nent2, nent3, nent4;
1796   Double_t int1, int2, int3, int4;
1797   nent1 = (Int_t)fLego->GetEntries();
1798   int1  = fLego->Integral();
1799   fC1->cd(1);
1800   if(int1) fLego->Draw("lego");
1801
1802   nent2 = (Int_t)fhLegoTracks->GetEntries();
1803   int2  = fhLegoTracks->Integral();
1804   fC1->cd(2);
1805   if(int2) fhLegoTracks->Draw("lego");
1806
1807   nent3 = (Int_t)fhLegoEMCAL->GetEntries();
1808   int3  = fhLegoEMCAL->Integral();
1809   fC1->cd(3);
1810   if(int3) fhLegoEMCAL->Draw("lego");
1811
1812   nent4 = (Int_t)fhLegoHadrCorr->GetEntries();
1813   int4  = fhLegoHadrCorr->Integral();
1814   fC1->cd(4);
1815   if(int4) fhLegoHadrCorr->Draw("lego");
1816
1817   // just for checking
1818   printf(" Integrals \n");
1819   printf("lego   %10.3f \ntrack  %10.3f \nhits   %10.3f \nHCorr   %10.3f\n--      %10.3f(must be 0)\n",
1820   int1, int2, int3, int4, int1 - (int2 + int3 - int4));
1821 }
1822
1823 const Char_t* AliEMCALJetFinder::GetFileNameForParameters(const char* dir)
1824 {
1825 // Get paramters from a file
1826   static TString tmp;
1827   if(strlen(dir)) tmp = dir;
1828   tmp += GetTitle();
1829   tmp += "_Par.txt";
1830   return tmp.Data();
1831 }
1832
1833 void AliEMCALJetFinder::RearrangeParticlesMemory(Int_t npart)
1834 {
1835         // See FillFromTracks() - npart must be positive
1836     if (fTrackList) delete[] fTrackList;
1837     if (fPtT)       delete[] fPtT;
1838     if (fEtaT)      delete[] fEtaT;
1839     if (fPhiT)      delete[] fPhiT;
1840     if (fPdgT)      delete[] fPdgT;
1841
1842     if(npart>0) {
1843        fTrackList = new Int_t  [npart];
1844        fPtT       = new Float_t[npart];
1845        fEtaT      = new Float_t[npart];
1846        fPhiT      = new Float_t[npart];
1847        fPdgT      = new Int_t[npart];
1848     } else {
1849        printf("AliEMCALJetFinder::RearrangeParticlesMemory : npart = %d\n", npart);
1850     }
1851 }
1852
1853 Bool_t AliEMCALJetFinder::IsThisPartonsOrDiQuark(Int_t pdg)
1854 {
1855 // Return quark info
1856   Int_t absPdg = TMath::Abs(pdg);
1857   if(absPdg<=6) return kTRUE; // quarks
1858   if(pdg == 21) return kTRUE; // gluon
1859   if(pdg == 92) return kTRUE; // string
1860
1861   // see p.51 of Pythia Manual
1862   // Not include diquarks with c and b quark - 4-mar-2002
1863   //                 ud_0,sd_0,su_0; dd_1,ud_1,uu_1;  sd_1,su_1,ss_1
1864   static Int_t diquark[9]={2101,3101,3201, 1103,2103,2203,  3103,3203,3303};
1865   for(Int_t i=0; i<9; i++) if(absPdg == diquark[i])  return kTRUE; // diquarks
1866
1867   return kFALSE;
1868 }
1869
1870 void AliEMCALJetFinder::FindChargedJet()
1871 {
1872 //
1873 // Find jet from charged particle information only
1874 //
1875
1876 //
1877 //  Look for seeds
1878 //
1879     Int_t njets = 0;
1880     Int_t part  = 0;
1881     Int_t nseed = 0;
1882
1883 //
1884 //
1885     ResetJets();
1886
1887 //
1888     for (part = 0; part < fNt; part++) {
1889         if (!fTrackList[part]) continue;
1890         if (fPtT[part] > fEtSeed) nseed++;
1891     }
1892     printf("\nThere are %d seeds (%d)\n", nseed, fNtS);
1893     Int_t* iSeeds = new Int_t[nseed];
1894     nseed = 0;
1895
1896     for (part = 0; part < fNt; part++) {
1897         if (!fTrackList[part]) continue;
1898         if (fPtT[part] > fEtSeed) iSeeds[nseed++] =  part;
1899     }
1900
1901 //
1902 // Loop over seeds
1903 //
1904     Int_t seed = 0;
1905     Float_t pt;
1906
1907     while(1){
1908 //
1909 // Find seed with highest pt
1910 //
1911         Float_t ptmax = -1.;
1912         Int_t   index = -1;
1913         Int_t   jndex = -1;
1914         for (seed = 0; seed < nseed; seed++) {
1915             if ((pt = fPtT[iSeeds[seed]]) > ptmax && iSeeds[seed] != -1) {
1916                 ptmax = pt;
1917                 index = seed;
1918             } // ptmax ?
1919         } // seeds
1920         if (ptmax < 0.) break;
1921         jndex = iSeeds[index];
1922 //
1923 // Remove from the list
1924         iSeeds[index] = -1;
1925         printf("\n Next Seed %d %f", jndex, ptmax);
1926 //
1927 // Find tracks in cone around seed
1928 //
1929         Float_t phiSeed = fPhiT[jndex];
1930         Float_t etaSeed = fEtaT[jndex];
1931         Float_t eT = 0.;
1932         Float_t pxJ = 0.;
1933         Float_t pyJ = 0.;
1934         Float_t pzJ = 0.;
1935
1936         for (part = 0; part < fNt; part++) {
1937             if (!fTrackList[part]) continue;
1938             Float_t deta = fEtaT[part] - etaSeed;
1939             Float_t dphi = fPhiT[part] - phiSeed;
1940             Float_t dR   = TMath::Sqrt(deta * deta + dphi * dphi);
1941             if (dR < fConeRadius) {
1942                 eT += fPtT[part];
1943                 Float_t theta = 2. * TMath::ATan(TMath::Exp(-fEtaT[part]));
1944                 Float_t px = fPtT[part] * TMath::Cos(fPhiT[part]);
1945                 Float_t py = fPtT[part] * TMath::Sin(fPhiT[part]);
1946                 Float_t pz = fPtT[part] / TMath::Tan(theta);
1947                 pxJ += px;
1948                 pyJ += py;
1949                 pzJ += pz;
1950                 //
1951                 // if seed, remove it
1952                 //
1953                 for (seed = 0; seed < nseed; seed++) {
1954                     if (part == iSeeds[seed]) iSeeds[seed] = -1;
1955                 } // seed ?
1956             } // < cone radius
1957         } // particle loop
1958
1959 //
1960 //      Estimate of jet direction
1961         Float_t phiJ   = TMath::ATan2(pyJ, pxJ);
1962         Float_t thetaJ = TMath::ATan2(TMath::Sqrt(pxJ * pxJ + pyJ * pyJ), pzJ);
1963         Float_t etaJ   = TMath::Log(TMath::Tan(thetaJ / 2.));
1964 //      Float_t ptJ    = TMath::Sqrt(pxJ * pxJ + pyJ * pyJ);
1965
1966 //
1967 //      Sum up all energy
1968 //
1969         Int_t iPhi0 = Int_t((phiJ - fPhiMin) / fDphi);
1970         Int_t iEta0 = Int_t((etaJ - fEtaMin) / fDeta);
1971         Int_t dIphi = Int_t(fConeRadius / fDphi);
1972         Int_t dIeta = Int_t(fConeRadius / fDeta);
1973         Int_t iPhi, iEta;
1974         Float_t sumE = 0;
1975         for (iPhi = iPhi0 -dIphi; iPhi < iPhi0 + dIphi; iPhi++) {
1976             for (iEta = iEta0 - dIeta; iEta < iEta0 + dIeta; iEta++) {
1977                 if (iPhi < 0 || iEta < 0) continue;
1978                 Float_t dPhi = fPhiMin + iPhi * fDphi;
1979                 Float_t dEta = fEtaMin + iEta * fDeta;
1980                 if (TMath::Sqrt(dPhi * dPhi + dEta * dEta) < fConeRadius) continue;
1981                 sumE += fLego->GetBinContent(iEta, iPhi);
1982             } // eta
1983         } // phi
1984 //
1985 //
1986 //
1987         fJetT[njets++] = new AliEMCALJet(sumE, phiJ, etaJ);
1988         FindTracksInJetCone();
1989         printf("\n Jet Energy %f %f %f %f %d\n", eT, sumE, fPtT[6], fPtT[7], njets);
1990         printf("\n Jet Phi %f %f %f \n", phiJ, fPhiT[6], fPhiT[7]);
1991         printf("\n Jet Eta %f %f %f \n", etaJ, fEtaT[6], fEtaT[7]);
1992     } // while(1)
1993     EMCALJETS.njet = njets;
1994     if (fWrite) WriteJets();
1995     fEvent++;
1996 }
1997 // 16-jan-2003 - just for convenience
1998 void AliEMCALJetFinder::Browse(TBrowser* b)
1999 {
2000 // Add to browser
2001    if(fHistsList)  b->Add((TObject*)fHistsList);
2002 }
2003
2004 Bool_t AliEMCALJetFinder::IsFolder() const
2005 {
2006 // Return folder status
2007   if(fHistsList) return kTRUE;
2008   else           return kFALSE;
2009 }
2010
2011 const Char_t* AliEMCALJetFinder::GetNameOfVariant()
2012 {
2013         // generate the literal string with info about jet finder
2014   Char_t name[200];
2015   sprintf(name, "jF_R%3.2fMinCell%4.1fPtCut%4.1fEtSeed%4.1fMinEt%4.1fBGSubtr%iSF%4.1f",
2016           fConeRadius,fMinCellEt,fPtCut,fEtSeed,fMinJetEt, fMode, fSamplingF);
2017   TString nt(name);
2018   nt.ReplaceAll(" ","");
2019   if(fBGFileName.Length()) {
2020      Int_t i1 = fBGFileName.Index("kBackground");
2021      Int_t i2 = fBGFileName.Index("/0000") - 1;
2022      if(i1>=0 && i2>=0) {
2023         TString bg(fBGFileName(i1,i2-i1+1));
2024         nt += bg;
2025      }
2026   }
2027   printf("<I> Name of variant %s \n", nt.Data());
2028   return nt.Data();
2029 }