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 /* $Id$ */
  17 //*-- Author: Andreas Morsch (CERN)
  18 #include <TClonesArray.h>
  19 #include <TTree.h>
  20 #include <TFile.h>
  21 #include <TH2.h>
  22 #include <TAxis.h>
  23 #include <TParticle.h>
  24 #include <TParticlePDG.h>
  25 #include "AliEMCALJetFinder.h"
  26 #include "AliEMCALGeometry.h"
  27 #include "AliEMCALHit.h"
  28 #include "Ecommon.h"
  29 #include "AliRun.h"
  30 #include "AliEMCAL.h"
  31 #include "AliHeader.h"
  32
  33 ClassImp(AliEMCALJetFinder)
  34
  35 //____________________________________________________________________________
  36 AliEMCALJetFinder::AliEMCALJetFinder()
  37 {
  38 // Default constructor
  39     fJets      = 0;
  40     fNjets     = 0;
  41     fLego      = 0;
  42     fTrackList = 0;
  43     fPtT       = 0;
  44     fEtaT      = 0;
  45     fPhiT      = 0;
  46 }
  47
  48 AliEMCALJetFinder::AliEMCALJetFinder(const char* name, const char *title)
  49     : TTask(name, title)
  50 {
  51 // Constructor
  52     fJets  = new TClonesArray("AliEMCALJet",10000);
  53     fNjets = 0;
  54     for (Int_t i=0; i<30000; i++)
  55     {
  56         fEtCell[i]  = 0.;
  57         fEtaCell[i] = 0.;
  58         fPhiCell[i] = 0.;
  59     }
  60     fLego = 0;
  61     fTrackList = 0;
  62     fTrackList = 0;
  63     fPtT       = 0;
  64     fEtaT      = 0;
  65     fPhiT      = 0;
  66 }
  67
  68
  69
  70
  71 //____________________________________________________________________________
  72 AliEMCALJetFinder::~AliEMCALJetFinder()
  73 {
  74 // Destructor
  75 //
  76     if (fJets){
  77         fJets->Delete();
  78         delete fJets;
  79     }
  80 }
  81
  82
  83
  84
  85 #ifndef WIN32
  86 # define jet_finder_ua1 jet_finder_ua1_
  87 # define hf1 hf1_
  88 # define type_of_call
  89
  90 #else
  91 # define jet_finder_ua1 J
  92 # define hf1 HF1
  93 # define type_of_call _stdcall
  94 #endif
  95
  96 extern "C" void type_of_call
  97 jet_finder_ua1(Int_t& ncell, Int_t& ncell_tot,
  98                Float_t  etc[30000],  Float_t etac[30000],
  99                Float_t  phic[30000],
 100                Float_t& min_move, Float_t& max_move, Int_t& mode,
 101                Float_t& prec_bg,  Int_t& ierror);
 102
 103 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
 104
 105
 106
 107
 108
 109 void AliEMCALJetFinder::Find(Int_t ncell, Int_t ncell_tot, Float_t etc[30000],
 110                              Float_t etac[30000], Float_t phic[30000],
 111                              Float_t min_move, Float_t max_move, Int_t mode,
 112                              Float_t prec_bg,  Int_t   ierror)
 113 {
 114 // Wrapper for fortran coded jet finder
 115 // Full argument list
 116     jet_finder_ua1(ncell, ncell_tot, etc, etac, phic,
 117                    min_move, max_move, mode, prec_bg, ierror);
 118     // Write result to output
 119     WriteJets();
 120 }
 121
 122 void AliEMCALJetFinder::Find()
 123 {
 124 // Wrapper for fortran coded jet finder using member data for
 125 // argument list
 126
 127     Float_t min_move = 0;
 128     Float_t max_move = 0;
 129     Int_t   mode     = 0;
 130     Float_t prec_bg  = 0.;
 131     Int_t   ierror   = 0;
 132
 133     jet_finder_ua1(fNcell, fNtot, fEtCell, fEtaCell, fPhiCell,
 134                    min_move, max_move, mode, prec_bg, ierror);
 135     // Write result to output
 136     Int_t njet = Njets();
 137     for (Int_t nj=0; nj<njet; nj++)
 138     {
 139         fJetT[nj] = new AliEMCALJet(JetEnergy(nj),
 140                                     JetPhiW(nj),
 141                                     JetEtaW(nj));
 142     }
 143     FindTracksInJetCone();
 144     WriteJets();
 145 }
 146
 147
 148 Int_t AliEMCALJetFinder::Njets()
 149 {
 150 // Get number of reconstructed jets
 151     return EMCALJETS.njet;
 152 }
 153
 154 Float_t AliEMCALJetFinder::JetEnergy(Int_t i)
 155 {
 156 // Get reconstructed jet energy
 157     return EMCALJETS.etj[i];
 158 }
 159
 160 Float_t AliEMCALJetFinder::JetPhiL(Int_t i)
 161 {
 162 // Get reconstructed jet phi from leading particle
 163     return EMCALJETS.phij[0][i];
 164 }
 165
 166 Float_t AliEMCALJetFinder::JetPhiW(Int_t i)
 167 {
 168 // Get reconstructed jet phi from weighting
 169     return EMCALJETS.phij[1][i];
 170 }
 171
 172 Float_t  AliEMCALJetFinder::JetEtaL(Int_t i)
 173 {
 174 // Get reconstructed jet eta from leading particles
 175     return EMCALJETS.etaj[0][i];
 176 }
 177
 178
 179 Float_t  AliEMCALJetFinder::JetEtaW(Int_t i)
 180 {
 181 // Get reconstructed jet eta from weighting
 182     return EMCALJETS.etaj[1][i];
 183 }
 184
 185 void AliEMCALJetFinder::SetCellSize(Float_t eta, Float_t phi)
 186 {
 187 // Set grid cell size
 188     EMCALCELLGEO.etaCellSize = eta;
 189     EMCALCELLGEO.phiCellSize = phi;
 190 }
 191
 192 void AliEMCALJetFinder::SetConeRadius(Float_t par)
 193 {
 194 // Set jet cone radius
 195     EMCALJETPARAM.coneRad = par;
 196     fConeRadius = par;
 197 }
 198
 199 void AliEMCALJetFinder::SetEtSeed(Float_t par)
 200 {
 201 // Set et cut for seeds
 202     EMCALJETPARAM.etSeed = par;
 203 }
 204
 205 void AliEMCALJetFinder::SetMinJetEt(Float_t par)
 206 {
 207 // Set minimum jet et
 208     EMCALJETPARAM.ejMin = par;
 209 }
 210
 211 void AliEMCALJetFinder::SetMinCellEt(Float_t par)
 212 {
 213 // Set et cut per cell
 214     EMCALJETPARAM.etMin = par;
 215 }
 216
 217 void AliEMCALJetFinder::SetPtCut(Float_t par)
 218 {
 219 // Set pT cut on charged tracks
 220     fPtCut = par;
 221 }
 222
 223
 224 void AliEMCALJetFinder::Test()
 225 {
 226 //
 227 // Test the finder call
 228 //
 229     const Int_t nmax = 30000;
 230     Int_t ncell      = 10;
 231     Int_t ncell_tot  = 100;
 232
 233     Float_t etc[nmax];
 234     Float_t etac[nmax];
 235     Float_t phic[nmax];
 236     Float_t min_move = 0;
 237     Float_t max_move = 0;
 238     Int_t   mode = 0;
 239     Float_t prec_bg = 0;
 240     Int_t   ierror = 0;
 241
 242
 243     Find(ncell, ncell_tot, etc, etac, phic,
 244          min_move, max_move, mode, prec_bg, ierror);
 245
 246 }
 247
 248 //
 249 //  I/O
 250 //
 251
 252 void AliEMCALJetFinder::AddJet(const AliEMCALJet& jet)
 253 {
 254     //
 255     // Add a jet
 256     //
 257     TClonesArray &lrawcl = *fJets;
 258     new(lrawcl[fNjets++]) AliEMCALJet(jet);
 259 }
 260
 261 void AliEMCALJetFinder::ResetJets()
 262 {
 263     //
 264     // Reset Jet List
 265     //
 266     fJets->Clear();
 267     fNjets = 0;
 268 }
 269
 270 void AliEMCALJetFinder::WriteJets()
 271 {
 272 //
 273 // Add all jets to the list
 274 //
 275     const Int_t kBufferSize = 4000;
 276     TTree *pK = gAlice->TreeK();
 277     const char* file = (pK->GetCurrentFile())->GetName();
 278 // I/O
 279     AliEMCAL* pEMCAL = (AliEMCAL* )gAlice->GetModule("EMCAL");
 280     printf("Make Branch - TreeR address %p %p\n",gAlice->TreeR(), pEMCAL);
 281     if (fJets && gAlice->TreeR()) {
 282         pEMCAL->MakeBranchInTree(gAlice->TreeR(),
 283                                  "Jets",
 284                                  &fJets,
 285                                  kBufferSize,
 286                                  file);
 287     }
 288     Int_t njet = Njets();
 289     for (Int_t nj = 0; nj < njet; nj++)
 290     {
 291         AddJet(*fJetT[nj]);
 292         delete fJetT[nj];
 293     }
 294
 295     Int_t nev = gAlice->GetHeader()->GetEvent();
 296     gAlice->TreeR()->Fill();
 297     char hname[30];
 298     sprintf(hname,"TreeR%d", nev);
 299     gAlice->TreeR()->Write(hname);
 300     gAlice->TreeR()->Reset();
 301     ResetJets();
 302 }
 303
 304 void AliEMCALJetFinder::BookLego()
 305 {
 306 //
 307 //  Book histo for discretisation
 308 //
 309 //
 310 //  Get geometry parameters from
 311     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 312     AliEMCALGeometry* geom =
 313         AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 314     fNbinEta = geom->GetNZ();
 315     fNbinPhi = geom->GetNPhi();
 316     const Float_t  phiMin  = geom->GetArm1PhiMin()*TMath::Pi()/180.;
 317     const Float_t  phiMax  = geom->GetArm1PhiMax()*TMath::Pi()/180.;
 318     fDphi   = (phiMax-phiMin)/fNbinEta;
 319     fDeta   = 1.4/fNbinEta;
 320     fNtot   = fNbinPhi*fNbinEta;
 321
 322 //
 323     fLego = new TH2F("legoH","eta-phi",
 324                            fNbinEta, -0.7,  0.7,
 325                            fNbinPhi, phiMin, phiMax);
 326 }
 327
 328 void AliEMCALJetFinder::DumpLego()
 329 {
 330 //
 331 // Dump lego histo into array
 332 //
 333     fNcell = 0;
 334     for (Int_t i = 1; i < fNbinEta; i++) {
 335         for (Int_t j = 1; j < fNbinPhi; j++) {
 336             Float_t e    = fLego->GetBinContent(i,j);
 337             TAxis* Xaxis = fLego->GetXaxis();
 338             TAxis* Yaxis = fLego->GetYaxis();
 339             Float_t eta  = Xaxis->GetBinCenter(i);
 340             Float_t phi  = Yaxis->GetBinCenter(j);
 341             fEtCell[fNcell]  = e;
 342             fEtaCell[fNcell] = eta;
 343             fPhiCell[fNcell] = phi;
 344             fNcell++;
 345         } // phi
 346     } // eta
 347     fNcell--;
 348 }
 349
 350 void AliEMCALJetFinder::ResetMap()
 351 {
 352 //
 353 // Reset eta-phi array
 354
 355     for (Int_t i=0; i<30000; i++)
 356     {
 357         fEtCell[i]  = 0.;
 358         fEtaCell[i] = 0.;
 359         fPhiCell[i] = 0.;
 360     }
 361 }
 362
 363
 364 void AliEMCALJetFinder::FillFromTracks(Int_t flag, Int_t ich)
 365 {
 366 //
 367 // Fill Cells with hit information
 368 //
 369 //
 370     ResetMap();
 371
 372     if (!fLego) BookLego();
 373 // Reset
 374     if (flag == 0) fLego->Reset();
 375 //
 376 // Access particle information
 377     Int_t npart = (gAlice->GetHeader())->GetNprimary();
 378 // Create track list
 379 //
 380 // 0: not selected
 381 // 1: selected for jet finding
 382 // 2: inside jet
 383 // ....
 384     if (fTrackList) delete[] fTrackList;
 385     if (fPtT)       delete[] fPtT;
 386     if (fEtaT)      delete[] fEtaT;
 387     if (fPhiT)      delete[] fPhiT;
 388
 389     fTrackList = new Int_t  [npart];
 390     fPtT       = new Float_t[npart];
 391     fEtaT      = new Float_t[npart];
 392     fPhiT      = new Float_t[npart];
 393     fNt        = npart;
 394
 395     for (Int_t part = 0; part < npart; part++) {
 396         TParticle *MPart = gAlice->Particle(part);
 397         Int_t mpart   = MPart->GetPdgCode();
 398         Int_t child1  = MPart->GetFirstDaughter();
 399         Float_t pT    = MPart->Pt();
 400         Float_t phi   = MPart->Phi();
 401         Float_t eta   = MPart->Eta();
 402 //      if (part == 6 || part == 7)
 403 //      {
 404 //          printf("\n Simulated Jet (pt, eta, phi): %d %f %f %f",
 405 //                 part-5, pT, eta, phi);
 406 //      }
 407
 408         fTrackList[part] = 0;
 409         fPtT[part]       = pT;
 410         fEtaT[part]      = eta;
 411         fPhiT[part]      = phi;
 412
 413         if (part < 2) continue;
 414         if (pT == 0 || pT < fPtCut) continue;
 415 // charged or neutral
 416         if (ich == 0) {
 417             TParticlePDG* pdgP = MPart->GetPDG();
 418             if (pdgP->Charge() == 0) continue;
 419         }
 420 // skip partons
 421         if (TMath::Abs(mpart) <= 6         ||
 422             mpart == 21                    ||
 423             mpart == 92) continue;
 424 // acceptance cut
 425         if (TMath::Abs(eta) > 0.7)         continue;
 426         if (phi*180./TMath::Pi() > 120.)   continue;
 427 // final state only
 428         if (child1 >= 0 && child1 < npart) continue;
 429 //      printf("\n sel:%5d %5d %5d %8.2f %8.2f %8.2f",
 430 //      part, mpart, child1, eta, phi, pT);
 431         fTrackList[part] = 1;
 432         fLego->Fill(eta, phi, pT);
 433     } // primary loop
 434     DumpLego();
 435 }
 436
 437 void AliEMCALJetFinder::FillFromHits(Int_t flag)
 438 {
 439 //
 440 // Fill Cells with track information
 441 //
 442 //
 443     ResetMap();
 444
 445     if (!fLego) BookLego();
 446     if (flag == 0) fLego->Reset();
 447
 448 //
 449 // Access hit information
 450     AliEMCAL* pEMCAL = (AliEMCAL*) gAlice->GetModule("EMCAL");
 451 //    AliEMCALGeometry* geom =
 452 //      AliEMCALGeometry::GetInstance(pEMCAL->GetTitle(), "");
 453
 454     TTree *treeH = gAlice->TreeH();
 455     Int_t ntracks = (Int_t) treeH->GetEntries();
 456 //
 457 //   Loop over tracks
 458 //
 459     Int_t nbytes = 0;
 460
 461
 462     for (Int_t track=0; track<ntracks;track++) {
 463         gAlice->ResetHits();
 464         nbytes += treeH->GetEvent(track);
 465 //
 466 //  Loop over hits
 467 //
 468         for(AliEMCALHit* mHit=(AliEMCALHit*) pEMCAL->FirstHit(-1);
 469             mHit;
 470             mHit=(AliEMCALHit*) pEMCAL->NextHit())
 471         {
 472             Float_t x      =    mHit->X();         // x-pos of hit
 473             Float_t y      =    mHit->Y();         // y-pos
 474             Float_t z      =    mHit->Z();         // z-pos
 475             Float_t eloss  =    mHit->GetEnergy(); // deposited energy
 476 //          Int_t   index  =    mHit->GetId();     // cell index
 477 //          Float_t eta, phi;
 478 //          geom->EtaPhiFromIndex(index,  eta, phi);
 479             Float_t r      = TMath::Sqrt(x*x+y*y);
 480             Float_t theta  = TMath::ATan2(r,z);
 481             Float_t eta    = -TMath::Log(TMath::Tan(theta/2.));
 482             Float_t phi    = TMath::ATan2(y,x);
 483             fLego->Fill(eta, phi, eloss);
 484 //          if (eloss > 1.) printf("\nx,y,z %f %f %f %f %f",
 485 //          r, z, eta, phi, eloss);
 486 //          printf("\n Max %f", fLego->GetMaximum());
 487         } // Hit Loop
 488     } // Track Loop
 489     DumpLego();
 490 }
 491
 492 void AliEMCALJetFinder::FindTracksInJetCone()
 493 {
 494 //
 495 //  Build list of tracks inside jet cone
 496 //
 497 //  Loop over jets
 498     Int_t njet = Njets();
 499     for (Int_t nj = 0; nj < njet; nj++)
 500     {
 501         Float_t etaj = JetEtaW(nj);
 502         Float_t phij = JetPhiW(nj);
 503         Int_t   nT   = 0;           // number of associated tracks
 504
 505 // Loop over particles
 506         for (Int_t part = 0; part < fNt; part++) {
 507             if (!fTrackList[part]) continue;
 508             Float_t phi      = fPhiT[part];
 509             Float_t eta      = fEtaT[part];
 510             Float_t dr       = TMath::Sqrt((etaj-eta)*(etaj-eta) +
 511                                            (phij-phi)*(phij-phi));
 512             if (dr < fConeRadius) {
 513                 fTrackList[part] = nj+2;
 514                 nT++;
 515             } // < ConeRadius ?
 516         } // particle loop
 517         Float_t* ptT  = new Float_t[nT];
 518         Float_t* etaT = new Float_t[nT];
 519         Float_t* phiT = new Float_t[nT];
 520         Int_t    iT   = 0;
 521         for (Int_t part = 0; part < fNt; part++) {
 522             if (fTrackList[part] == nj+2) {
 523                 ptT [iT] = fPtT [part];
 524                 etaT[iT] = fEtaT[part];
 525                 phiT[iT] = fPhiT[part];
 526                 iT++;
 527             } // particle associated
 528         } // particle loop
 529         fJetT[nj]->SetTrackList(nT, ptT, etaT, phiT);
 530         delete[] ptT;
 531         delete[] etaT;
 532         delete[] phiT;
 533     } // jet loop loop
 534 }
 535
 536
 537
 538 void hf1(Int_t& id, Float_t& x, Float_t& wgt)
 539 {
 540 // dummy for hbook calls
 541     ;
 542 }
 543
 544