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