EMCAL/AliEMCALJetFinder.cxx

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