JETAN/AliUA1JetFinder.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id */
  17
  18 //---------------------------------------------------------------------
  19 // UA1 Cone Algorithm Jet finder for jet studies
  20 // manages the search for jets using charged particle momentum and
  21 // neutral cell energy information
  22 // Authors: Rafael.Diaz.Valdes@cern.ch
  23 //          magali.estienne@subatech.in2p3.fr
  24 //          alexandre.shabetai@cern.ch
  25 // ** 2011
  26 // Modified accordingly to reader/finder splitting and new handling of neutral information
  27 // Versions V1 and V2 merged
  28 //---------------------------------------------------------------------
  29
  30 #include <TH2F.h>
  31 #include <TMath.h>
  32
  33 #include "AliUA1JetFinder.h"
  34 #include "AliUA1JetHeaderV1.h"
  35 #include "AliJetCalTrk.h"
  36 #include "AliJetBkg.h"
  37 #include "AliAODJetEventBackground.h"
  38 #include "AliAODJet.h"
  39
  40 ClassImp(AliUA1JetFinder)
  41
  42 ////////////////////////////////////////////////////////////////////////
  43
  44 AliUA1JetFinder::AliUA1JetFinder():
  45   AliJetFinder(),
  46   fLego(0),
  47   fJetBkg(new AliJetBkg())
  48 {
  49   // Default constructor
  50 }
  51
  52 //-----------------------------------------------------------------------
  53 AliUA1JetFinder::~AliUA1JetFinder()
  54 {
  55   // Destructor
  56   delete fLego;
  57   delete fJetBkg;
  58
  59 }
  60
  61 //-----------------------------------------------------------------------
  62 void AliUA1JetFinder::FindJets()
  63 {
  64   //  Used to find jets using charged particle momentum information
  65   //  & neutral energy from calo cells
  66   //
  67   //  1) Fill cell map array
  68   //  2) calculate total energy and fluctuation level
  69   //  3) Run algorithm
  70   //     3.1) look centroides in cell map
  71   //     3.2) calculate total energy in cones
  72   //     3.3) flag as a possible jet
  73   //     3.4) reorder cones by energy
  74   //  4) subtract backg in accepted jets
  75   //  5) fill AliJet list
  76
  77   // transform input to pt,eta,phi plus lego
  78
  79   AliUA1JetHeader* header  = (AliUA1JetHeader*) fHeader;
  80   Int_t              nIn = fCalTrkEvent->GetNCalTrkTracks();
  81   fDebug   = fHeader->GetDebug();
  82
  83   if (nIn <= 0) return;
  84   fJetBkg->SetHeader(fHeader);
  85   fJetBkg->SetCalTrkEvent(GetCalTrkEvent());
  86   fJetBkg->SetDebug(fDebug);
  87   // local arrays for input
  88   // ToDo: check memory fragmentation, maybe better to
  89   // define them globally and resize as needed
  90   // Fragmentation should be worse for low mult...
  91   Float_t* ptT      = new Float_t[nIn];
  92   Float_t* etaT     = new Float_t[nIn];
  93   Float_t* phiT     = new Float_t[nIn];
  94   Int_t*   injet    = new Int_t[nIn];
  95   Int_t*   injetOk  = new Int_t[nIn];
  96
  97   memset(ptT,0,sizeof(Float_t)*nIn);
  98   memset(etaT,0,sizeof(Float_t)*nIn);
  99   memset(phiT,0,sizeof(Float_t)*nIn);
 100   memset(injet,0,sizeof(Int_t)*nIn);
 101   memset(injetOk,-1,sizeof(Int_t)*nIn);
 102
 103   // load input vectors and calculate total energy in array
 104
 105   // total energy in array
 106   Float_t etbgTotal = 0.;
 107   Float_t npart = 0.;
 108   Float_t etbg2 = 0.;
 109
 110   for (Int_t i = 0; i < fCalTrkEvent->GetNCalTrkTracks(); i++){
 111     ptT[i]  = fCalTrkEvent->GetCalTrkTrack(i)->GetPt();
 112     etaT[i] = fCalTrkEvent->GetCalTrkTrack(i)->GetEta();
 113     phiT[i] = ((fCalTrkEvent->GetCalTrkTrack(i)->GetPhi() < 0) ? (fCalTrkEvent->GetCalTrkTrack(i)->GetPhi()) + 2 * TMath::Pi() : fCalTrkEvent->GetCalTrkTrack(i)->GetPhi());
 114     //fCalTrkEvent->GetCalTrkTrack(i)->Print(Form("%d",i));
 115     if (fCalTrkEvent->GetCalTrkTrack(i)->GetCutFlag() != 1) continue;
 116     fLego->Fill(etaT[i], phiT[i], ptT[i]);
 117     npart += 1;
 118     etbgTotal+= ptT[i];
 119     etbg2 += ptT[i]*ptT[i];
 120   }
 121
 122   // calculate total energy and fluctuation in map
 123   Double_t meanpt = 0.;
 124   Double_t ptRMS = 0.;
 125   if(npart>0){
 126     meanpt = etbgTotal/npart;
 127     etbg2 = etbg2/npart;
 128     if(etbg2>(meanpt*meanpt)){// prenent NAN, should only happen due to numerical instabilities
 129       ptRMS = TMath::Sqrt(etbg2-meanpt*meanpt);
 130     }
 131   }
 132   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
 133
 134   // arrays to hold jets
 135   Float_t etaJet[kMaxJets];
 136   Float_t phiJet[kMaxJets];
 137   Float_t etJet[kMaxJets];
 138   Float_t etsigJet[kMaxJets]; //signal et in jet
 139   Float_t etallJet[kMaxJets];  // total et in jet (tmp variable)
 140   Int_t   ncellsJet[kMaxJets];
 141   Int_t   multJetT[kMaxJets];
 142   Int_t   multJetC[kMaxJets];
 143   Int_t   multJet[kMaxJets];
 144   Float_t *areaJet = new Float_t[kMaxJets];
 145   // Used for jet reordering at the end of the jet finding procedure
 146   Float_t etaJetOk[kMaxJets];
 147   Float_t phiJetOk[kMaxJets];
 148   Float_t etJetOk[kMaxJets];
 149   Float_t etsigJetOk[kMaxJets]; //signal et in jet
 150   Float_t etallJetOk[kMaxJets];  // total et in jet (tmp variable)
 151   Int_t   ncellsJetOk[kMaxJets];
 152   Int_t   multJetOk[kMaxJets];
 153   Float_t *areaJetOk = new Float_t[kMaxJets];
 154   Int_t   nJets; // to hold number of jets found by algorithm
 155   Int_t   nj;    // number of jets accepted
 156   Float_t prec  = header->GetPrecBg();
 157   Float_t bgprec = 1;
 158
 159   while(bgprec > prec){
 160     //reset jet arrays in memory
 161     memset(etaJet,0,sizeof(Float_t)*kMaxJets);
 162     memset(phiJet,0,sizeof(Float_t)*kMaxJets);
 163     memset(etJet,0,sizeof(Float_t)*kMaxJets);
 164     memset(etallJet,0,sizeof(Float_t)*kMaxJets);
 165     memset(etsigJet,0,sizeof(Float_t)*kMaxJets);
 166     memset(ncellsJet,0,sizeof(Int_t)*kMaxJets);
 167     memset(multJetT,0,sizeof(Int_t)*kMaxJets);
 168     memset(multJetC,0,sizeof(Int_t)*kMaxJets);
 169     memset(multJet,0,sizeof(Int_t)*kMaxJets);
 170     memset(areaJet,0,sizeof(Float_t)*kMaxJets);
 171     memset(etaJetOk,0,sizeof(Float_t)*kMaxJets);
 172     memset(phiJetOk,0,sizeof(Float_t)*kMaxJets);
 173     memset(etJetOk,0,sizeof(Float_t)*kMaxJets);
 174     memset(etallJetOk,0,sizeof(Float_t)*kMaxJets);
 175     memset(etsigJetOk,0,sizeof(Float_t)*kMaxJets);
 176     memset(ncellsJetOk,0,sizeof(Int_t)*kMaxJets);
 177     memset(multJetOk,0,sizeof(Int_t)*kMaxJets);
 178     memset(areaJetOk,0,sizeof(Float_t)*kMaxJets);
 179     nJets = 0;
 180     nj = 0;
 181     // reset particles-jet array in memory
 182     memset(injet,-1,sizeof(Int_t)*nIn);
 183     //run cone algorithm finder
 184     RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
 185     //run background subtraction
 186     if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
 187       nj = header->GetNAcceptJets();
 188     else
 189       nj = nJets;
 190
 191     //subtract background
 192     Float_t etbgTotalN = 0.0; //new background
 193     Float_t sigmaN = 0.0; //new background
 194     if(header->GetBackgMode() == 1) {// standard
 195       fJetBkg->SubtractBackg(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);}
 196     if(header->GetBackgMode() == 2) //cone
 197       fJetBkg->SubtractBackgCone(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
 198     if(header->GetBackgMode() == 3) //ratio
 199       fJetBkg->SubtractBackgRatio(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
 200     if(header->GetBackgMode() == 4) //statistic
 201       fJetBkg->SubtractBackgStat(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
 202
 203     //calc precision
 204     if(etbgTotalN != 0.0)
 205       bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
 206     else
 207       bgprec = 0;
 208     etbgTotal = etbgTotalN; // update with new background estimation
 209
 210   } //end while
 211
 212   // add tracks to the jet if it wasn't yet done
 213   if (header->GetBackgMode() == 0){
 214     Float_t rc= header->GetRadius();
 215     for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 216       for(Int_t ijet=0; ijet<nJets; ijet++){
 217         Float_t deta = etaT[jpart] - etaJet[ijet];
 218         Float_t dphi = phiT[jpart] - phiJet[ijet];
 219         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 220         if (dphi >  TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 221         Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 222         if(dr <= rc){ // particles inside this cone
 223           injet[jpart] = ijet;
 224           break;
 225         }
 226       }// end jets loop
 227     } //end particle loop
 228   }
 229
 230   // add jets to list
 231   if (fDebug>1) printf("Found %d jets \n", nj);
 232
 233   // Reorder jets by et in cone
 234   // Sort jets by energy
 235   Int_t idx[kMaxJets];
 236   TMath::Sort(nJets, etJet, idx);
 237   for(Int_t p = 0; p < nJets; p++)
 238     {
 239       etaJetOk[p]    = etaJet[idx[p]];
 240       phiJetOk[p]    = phiJet[idx[p]];
 241       etJetOk[p]     = etJet[idx[p]];
 242       etallJetOk[p]  = etJet[idx[p]];
 243       etsigJetOk[p]  = etsigJet[idx[p]];
 244       ncellsJetOk[p] = ncellsJet[idx[p]];
 245       multJetOk[p]   = multJet[idx[p]];
 246       areaJetOk[p]   = areaJet[idx[p]];
 247     }
 248
 249   //////////////////////////
 250
 251   Int_t nTracks = fCalTrkEvent->GetNCalTrkTracks();
 252
 253   for(Int_t kj=0; kj<nj; kj++)
 254     {
 255       if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
 256           (etaJetOk[kj] < (header->GetJetEtaMin())) ||
 257           (phiJetOk[kj] > (header->GetJetPhiMax())) ||
 258           (phiJetOk[kj] < (header->GetJetPhiMin())) ||
 259           (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
 260       Float_t px=-999, py=-999 ,pz=-999 ,en=-999; // convert to 4-vector
 261       px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
 262       py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
 263       pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
 264       en = TMath::Sqrt(px * px + py * py + pz * pz);
 265       AliAODJet jet(px, py, pz, en);
 266
 267       // Calc jet area if it wasn't yet done
 268       if (header->GetBackgMode() == 0){
 269         // calculate the area of the jet
 270         Float_t rc= header->GetRadius();
 271         areaJetOk[kj] = fJetBkg->CalcJetAreaEtaCut(rc,etaJetOk[kj]);
 272       }
 273
 274       jet.SetEffArea(areaJetOk[kj],0.,0.,0.);
 275       jet.SetBgEnergy(etbgTotal,0.);
 276       if (fDebug>1) jet.Print(Form("%d",kj));
 277
 278       for(Int_t jpart = 0; jpart < nTracks; jpart++) { // loop for all particles in array
 279         // Track to jet reordering
 280         if(injet[jpart] == idx[kj]){
 281           injetOk[jpart] = kj;
 282         }
 283         // Check if the particle belongs to the jet and add the ref
 284         if(injetOk[jpart] == kj && fCalTrkEvent->GetCalTrkTrack(jpart)->GetCutFlag() == 1) {
 285           jet.AddTrack(fCalTrkEvent->GetCalTrkTrack(jpart)->GetTrackObject());
 286         }
 287       }
 288
 289       AddJet(jet);
 290
 291     }
 292
 293   //delete
 294   delete[] ptT;
 295   delete[] etaT;
 296   delete[] phiT;
 297   delete[] injet;
 298   delete[] injetOk;
 299   delete[] areaJet;
 300   delete[] areaJetOk;
 301
 302 }
 303
 304 //-----------------------------------------------------------------------
 305 void AliUA1JetFinder::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
 306                                   Float_t* const etJet,Float_t* const etaJet, Float_t* const phiJet,
 307                                   Float_t* const etallJet, Int_t* const ncellsJet)
 308 {
 309   // Dump lego
 310   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
 311   const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory
 312
 313   const Int_t nBinEta = header->GetLegoNbinEta();
 314   const Int_t nBinPhi = header->GetLegoNbinPhi();
 315   if((nBinPhi*nBinEta)>nBinsMax){
 316     AliError("Too many bins of the ETA-PHI histogram");
 317   }
 318
 319   Float_t etCell[nBinsMax] = {0.};   //! Cell Energy
 320   Float_t etaCell[nBinsMax] = {0.};  //! Cell eta
 321   Float_t phiCell[nBinsMax] = {0.};  //! Cell phi
 322   Short_t flagCell[nBinsMax] = {0}; //! Cell flag
 323
 324   Int_t nCell = 0;
 325   TAxis* xaxis = fLego->GetXaxis();
 326   TAxis* yaxis = fLego->GetYaxis();
 327   Float_t e = 0.0;
 328   for (Int_t i = 1; i <= nBinEta; i++) {
 329     for (Int_t j = 1; j <= nBinPhi; j++) {
 330       e = fLego->GetBinContent(i,j);
 331       if (e < 0.0) continue; // don't include this cells
 332       Float_t eta  = xaxis->GetBinCenter(i);
 333       Float_t phi  = yaxis->GetBinCenter(j);
 334       etCell[nCell]  = e;
 335       etaCell[nCell] = eta;
 336       phiCell[nCell] = phi;
 337       flagCell[nCell] = 0; //default
 338       nCell++;
 339     }
 340   }
 341   // Parameters from header
 342   Float_t minmove = header->GetMinMove();
 343   Float_t maxmove = header->GetMaxMove();
 344   Float_t rc      = header->GetRadius();
 345   Float_t etseed  = header->GetEtSeed();
 346   // Tmp array of jets form algoritm
 347   Float_t etaAlgoJet[kMaxJets] = {0.0};
 348   Float_t phiAlgoJet[kMaxJets] = {0.0};
 349   Float_t etAlgoJet[kMaxJets] = {0.0};
 350   Int_t   ncellsAlgoJet[kMaxJets] = {0};
 351
 352   // Run algorithm//
 353
 354   // Sort cells by et
 355   Int_t index[nBinsMax];
 356   TMath::Sort(nCell, etCell, index);
 357   // variable used in centroide loop
 358   Float_t eta   = 0.0;
 359   Float_t phi   = 0.0;
 360   Float_t eta0  = 0.0;
 361   Float_t phi0  = 0.0;
 362   Float_t etab  = 0.0;
 363   Float_t phib  = 0.0;
 364   Float_t etas  = 0.0;
 365   Float_t phis  = 0.0;
 366   Float_t ets   = 0.0;
 367   Float_t deta  = 0.0;
 368   Float_t dphi  = 0.0;
 369   Float_t dr    = 0.0;
 370   Float_t etsb  = 0.0;
 371   Float_t etasb = 0.0;
 372   Float_t phisb = 0.0;
 373   Float_t dphib = 0.0;
 374
 375   for(Int_t icell = 0; icell < nCell; icell++)
 376     {
 377       Int_t jcell = index[icell];
 378       if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
 379       if(flagCell[jcell] != 0) continue; // if cell was used before
 380
 381       eta  = etaCell[jcell];
 382       phi  = phiCell[jcell];
 383       eta0 = eta;
 384       phi0 = phi;
 385       etab = eta;
 386       phib = phi;
 387       ets  = etCell[jcell];
 388       etas = 0.0;
 389       phis = 0.0;
 390       etsb = ets;
 391       etasb = 0.0;
 392       phisb = 0.0;
 393       for(Int_t kcell =0; kcell < nCell; kcell++)
 394         {
 395           Int_t lcell = index[kcell];
 396           if(lcell == jcell) continue; // cell itself
 397           if(flagCell[lcell] != 0) continue; // cell used before
 398           if(etCell[lcell] > etCell[jcell]) continue; // can this happen
 399           //calculate dr
 400           deta = etaCell[lcell] - eta;
 401           dphi = TMath::Abs(phiCell[lcell] - phi);
 402           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 403           dr = TMath::Sqrt(deta * deta + dphi * dphi);
 404           if(dr <= rc)
 405             {
 406               // calculate offset from initiate cell
 407               deta = etaCell[lcell] - eta0;
 408               dphi = phiCell[lcell] - phi0;
 409               if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
 410               if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
 411               etas = etas + etCell[lcell]*deta;
 412               phis = phis + etCell[lcell]*dphi;
 413               ets = ets + etCell[lcell];
 414               //new weighted eta and phi including this cell
 415               eta = eta0 + etas/ets;
 416               phi = phi0 + phis/ets;
 417               // if cone does not move much, just go to next step
 418               dphib = TMath::Abs(phi - phib);
 419               if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
 420               dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
 421               if(dr <= minmove) break;
 422               // cone should not move more than max_mov
 423               dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
 424               if(dr > maxmove){
 425                 eta = etab;
 426                 phi = phib;
 427                 ets = etsb;
 428                 etas = etasb;
 429                 phis = phisb;
 430               } else { // store this loop information
 431                 etab=eta;
 432                 phib=phi;
 433                 etsb = ets;
 434                 etasb = etas;
 435                 phisb = phis;
 436               }
 437             } // inside cone
 438         }//end of cells loop looking centroide
 439
 440       // Avoid cones overloap (to be implemented in the future)
 441
 442       // Flag cells in Rc, estimate total energy in cone
 443       Float_t etCone   = 0.0;
 444       Int_t   nCellIn  = 0;
 445       rc = header->GetRadius();
 446
 447       for(Int_t ncell =0; ncell < nCell; ncell++)
 448         {
 449           if(flagCell[ncell] != 0) continue; // cell used before
 450           //calculate dr
 451           deta = etaCell[ncell] - eta;
 452           dphi = phiCell[ncell] - phi;
 453           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 454           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 455           dr = TMath::Sqrt(deta * deta + dphi * dphi);
 456           if(dr <= rc){  // cell in cone
 457             flagCell[ncell] = -1;
 458             etCone+=etCell[ncell];
 459             nCellIn++;
 460           }
 461         }
 462
 463       // Select jets with et > background
 464       // estimate max fluctuation of background in cone
 465       Double_t ncellin = (Double_t)nCellIn;
 466       Double_t ntcell  = (Double_t)nCell;
 467       Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
 468       // min cone et
 469       Double_t etcmin = etCone ;  // could be used etCone - etmin !!
 470       //decisions !! etbmax < etcmin
 471
 472       for(Int_t mcell =0; mcell < nCell; mcell++){
 473         if(flagCell[mcell] == -1){
 474           if(etbmax < etcmin)
 475             flagCell[mcell] = 1; //flag cell as used
 476           else
 477             flagCell[mcell] = 0; // leave it free
 478         }
 479       }
 480       //store tmp jet info !!!
 481       if(etbmax < etcmin) {
 482         if(nJets<kMaxJets){
 483           etaAlgoJet[nJets]    = eta;
 484           phiAlgoJet[nJets]    = phi;
 485           etAlgoJet[nJets]     = etCone;
 486           ncellsAlgoJet[nJets] = nCellIn;
 487           nJets++;
 488         }
 489         else{
 490           AliError(Form("Too many jets (> %d) found by UA1JetFinder, adapt your cuts",kMaxJets));
 491           break;
 492         }
 493       }
 494     } // end of cells loop
 495
 496   //reorder jets by et in cone
 497   //sort jets by energy
 498   for(Int_t p = 0; p < nJets; p++)
 499     {
 500       etaJet[p]    = etaAlgoJet[p];
 501       phiJet[p]    = phiAlgoJet[p];
 502       etJet[p]     = etAlgoJet[p];
 503       etallJet[p]  = etAlgoJet[p];
 504       ncellsJet[p] = ncellsAlgoJet[p];
 505     }
 506
 507 }
 508
 509 //-----------------------------------------------------------------------
 510 void AliUA1JetFinder::Reset()
 511 {
 512   fLego->Reset();
 513   AliJetFinder::Reset();
 514
 515 }
 516
 517 //-----------------------------------------------------------------------
 518 void AliUA1JetFinder::WriteJHeaderToFile() const
 519 {
 520   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
 521   header->Write();
 522
 523 }
 524
 525 //-----------------------------------------------------------------------
 526 void AliUA1JetFinder::Init()
 527 {
 528
 529   // initializes some variables
 530   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
 531   // book lego
 532   fLego = new TH2F("legoH","eta-phi",
 533                    header->GetLegoNbinEta(), header->GetLegoEtaMin(),
 534                    header->GetLegoEtaMax(),  header->GetLegoNbinPhi(),
 535                    header->GetLegoPhiMin(),  header->GetLegoPhiMax());
 536
 537 }
 538