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