JETAN/AliUA1JetFinderV2.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
  20 // manages the search for jets
  21 // Author: Rafael.Diaz.Valdes@cern.ch
  22 // (version in c++)
  23 // Modified to include neutral particles (magali.estienne@ires.in2p3.fr)
  24 //---------------------------------------------------------------------
  25
  26 #include <Riostream.h>
  27
  28 #include <TArrayF.h>
  29 #include <TClonesArray.h>
  30 #include <TFile.h>
  31 #include <TH1F.h>
  32 #include <TH2F.h>
  33 #include <TLorentzVector.h>
  34 #include <TMath.h>
  35 #include <TRefArray.h>
  36
  37 #include "AliUA1JetFinderV2.h"
  38 #include "AliUA1JetHeaderV1.h"
  39 #include "AliJetUnitArray.h"
  40 #include "AliJetReaderHeader.h"
  41 #include "AliJetReader.h"
  42 #include "AliJet.h"
  43 #include "AliAODJet.h"
  44
  45
  46 ClassImp(AliUA1JetFinderV2)
  47
  48 ////////////////////////////////////////////////////////////////////////
  49
  50 AliUA1JetFinderV2::AliUA1JetFinderV2() :
  51     AliJetFinder(),
  52     fLego(0),
  53     fDebug(0),
  54     fOpt(0)
  55 {
  56   // Constructor
  57 }
  58
  59 ////////////////////////////////////////////////////////////////////////
  60
  61 AliUA1JetFinderV2::~AliUA1JetFinderV2()
  62
  63 {
  64   // destructor
  65 }
  66
  67 ////////////////////////////////////////////////////////////////////////
  68
  69
  70 void AliUA1JetFinderV2::FindJets()
  71
  72 {
  73   //1) Fill cell map array
  74   //2) calculate total energy and fluctuation level
  75   //3) Run algorithm
  76   //   3.1) look centroides in cell map
  77   //   3.2) calculate total energy in cones
  78   //   3.3) flag as a possible jet
  79   //   3.4) reorder cones by energy
  80   //4) subtract backg in accepted jets
  81   //5) fill AliJet list
  82
  83   // transform input to pt,eta,phi plus lego
  84
  85   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
  86   TClonesArray*      fUnit = fReader->GetUnitArray();
  87   Int_t              nCandidate = fReader->GetNumCandidate();
  88   Int_t              nIn = fUnit->GetEntries();
  89
  90   if (nIn == 0) return;
  91
  92   // local arrays for input No Cuts
  93   // Both pt < ptMin and pt > ptMin
  94   Float_t* enT       = new Float_t[nCandidate];
  95   Float_t* ptT       = new Float_t[nCandidate];
  96   Float_t* etaT      = new Float_t[nCandidate];
  97   Float_t* phiT      = new Float_t[nCandidate];
  98   Float_t* detaT     = new Float_t[nCandidate];
  99   Float_t* dphiT     = new Float_t[nCandidate];
 100   Float_t* cFlagT    = new Float_t[nCandidate];
 101   Float_t* cClusterT = new Float_t[nCandidate];
 102   Float_t* idT       = new Float_t[nCandidate];
 103   Int_t loop1 = 0;
 104   Int_t* injet = new Int_t[nCandidate];
 105   Int_t* sflag = new Int_t[nCandidate];
 106
 107
 108   //total energy in array
 109   Float_t  etbgTotal = 0.0;
 110   TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
 111
 112   // Input cell info
 113   Float_t *etCell = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
 114   Float_t *etaCell = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
 115   Float_t *phiCell = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
 116   Int_t   *flagCell = new Int_t[nIn];   //! Cell phi - Extracted from UnitArray
 117
 118   // Information extracted from fUnitArray
 119   // load input vectors and calculate total energy in array
 120   for(Int_t i=0; i<nIn; i++)
 121     {
 122       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
 123       if(uArray->GetUnitCutFlag()==1){
 124         etCell[i] = uArray->GetUnitEnergy();
 125         if (etCell[i] > 0.0) etCell[i] -= header->GetMinCellEt();
 126         if (etCell[i] < 0.0) etCell[i] = 0.;
 127         etaCell[i] = uArray->GetUnitEta();
 128         phiCell[i] = uArray->GetUnitPhi();
 129         flagCell[i] = 0; // default
 130       }
 131       else {
 132         etCell[i]  = 0.;
 133         etaCell[i] = uArray->GetUnitEta();
 134         phiCell[i] = uArray->GetUnitPhi();
 135         flagCell[i] = 0;
 136       }
 137
 138       if(uArray->GetUnitEnergy()>0.){
 139         ptT[loop1]   = uArray->GetUnitEnergy();
 140         enT[loop1]   = uArray->GetUnitEnergy();
 141         etaT[loop1]  = uArray->GetUnitEta();
 142         phiT[loop1]  = uArray->GetUnitPhi();
 143         detaT[loop1] = uArray->GetUnitDeta();
 144         dphiT[loop1] = uArray->GetUnitDphi();
 145         cFlagT[loop1]= uArray->GetUnitCutFlag();
 146         idT[loop1]   = uArray->GetUnitID();
 147         if(cFlagT[loop1] == 1) {
 148           hPtTotal->Fill(ptT[loop1]);
 149           //    fLego->Fill(etaT[i], phiT[i], ptT[i]);
 150           etbgTotal+= ptT[loop1];
 151         }
 152         loop1++;
 153       }
 154     }
 155
 156   //  fJets->SetNinput(nIn);
 157   fJets->SetNinput(nCandidate);
 158
 159   // calculate total energy and fluctuation in map
 160   Double_t meanpt = hPtTotal->GetMean();
 161   Double_t ptRMS = hPtTotal->GetRMS();
 162   Double_t npart = hPtTotal->GetEntries();
 163   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
 164
 165   // arrays to hold jets
 166   Float_t* etaJet = new Float_t[30];
 167   Float_t* phiJet = new Float_t[30];
 168   Float_t* etJet  = new Float_t[30];
 169   Float_t* etsigJet  = new Float_t[30]; //signal et in jet
 170   Float_t* etallJet = new Float_t[30];  // total et in jet (tmp variable)
 171   Int_t* ncellsJet = new Int_t[30];
 172   Int_t* multJet  = new Int_t[30];
 173   Int_t nJets; // to hold number of jets found by algorithm
 174   Int_t nj;    // number of jets accepted
 175   Float_t prec  = header->GetPrecBg();
 176   Float_t bgprec = 1;
 177   while(bgprec > prec){
 178      //reset jet arrays in memory
 179      memset(etaJet,0,sizeof(Float_t)*30);
 180      memset(phiJet,0,sizeof(Float_t)*30);
 181      memset(etJet,0,sizeof(Float_t)*30);
 182      memset(etallJet,0,sizeof(Float_t)*30);
 183      memset(etsigJet,0,sizeof(Float_t)*30);
 184      memset(ncellsJet,0,sizeof(Int_t)*30);
 185      memset(multJet,0,sizeof(Int_t)*30);
 186      nJets = 0;
 187      nj = 0;
 188      // reset particles-jet array in memory
 189      memset(injet,-1,sizeof(Int_t)*nCandidate);
 190      //run cone algorithm finder
 191      RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
 192      //run background subtraction
 193      if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
 194        nj = header->GetNAcceptJets();
 195      else
 196        nj = nJets;
 197      //subtract background
 198      Float_t etbgTotalN = 0.0; //new background
 199      if(header->GetBackgMode() == 1) // standar
 200         SubtractBackg(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 201      if(header->GetBackgMode() == 2) //cone
 202         SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 203      if(header->GetBackgMode() == 3) //ratio
 204         SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 205      if(header->GetBackgMode() == 4) //statistic
 206         SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 207      //calc precision
 208      if(etbgTotalN != 0.0)
 209         bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
 210      else
 211         bgprec = 0;
 212      etbgTotal = etbgTotalN; // update with new background estimation
 213   } //end while
 214
 215   // add jets to list
 216   Int_t* idxjets = new Int_t[nj];
 217   Int_t nselectj = 0;
 218   printf("Found %d jets \n", nj);
 219
 220   for(Int_t kj=0; kj<nj; kj++){
 221      if ((etaJet[kj] > (header->GetJetEtaMax())) ||
 222           (etaJet[kj] < (header->GetJetEtaMin())) ||
 223           (etJet[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
 224       Float_t px, py,pz,en; // convert to 4-vector
 225       px = etJet[kj] * TMath::Cos(phiJet[kj]);
 226       py = etJet[kj] * TMath::Sin(phiJet[kj]);
 227       pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj])));
 228       en = TMath::Sqrt(px * px + py * py + pz * pz);
 229       fJets->AddJet(px, py, pz, en);
 230       AliAODJet jet(px, py, pz, en);
 231       jet.Print("");
 232
 233       AddJet(jet);
 234
 235       idxjets[nselectj] = kj;
 236       nselectj++;
 237   }
 238   //add signal percentage and total signal  in AliJets for analysis tool
 239   Float_t* percentage  = new Float_t[nselectj];
 240   Int_t* ncells      = new Int_t[nselectj];
 241   Int_t* mult        = new Int_t[nselectj];
 242   for(Int_t i = 0; i< nselectj; i++){
 243      percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]];
 244      ncells[i] = ncellsJet[idxjets[i]];
 245      mult[i] = multJet[idxjets[i]];
 246   }
 247    //add particle-injet relationship ///
 248   //   for(Int_t bj = 0; bj < nIn; bj++){
 249    for(Int_t bj = 0; bj < nCandidate; bj++){
 250        if(injet[bj] == -1) continue; //background particle
 251        Int_t bflag = 0;
 252        for(Int_t ci = 0; ci< nselectj; ci++){
 253            if(injet[bj] == idxjets[ci]){
 254               injet[bj]= ci;
 255               bflag++;
 256               break;
 257            }
 258        }
 259        if(bflag == 0) injet[bj] = -1; // set as background particle
 260    }
 261   fJets->SetNCells(ncells);
 262   fJets->SetPtFromSignal(percentage);
 263   fJets->SetMultiplicities(mult);
 264   fJets->SetInJet(injet);
 265   fJets->SetEtaIn(etaT);
 266   fJets->SetPhiIn(phiT);
 267   fJets->SetPtIn(ptT);
 268   fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi()));
 269
 270
 271   //delete
 272   delete enT;
 273   delete ptT;
 274   delete etaT;
 275   delete phiT;
 276   delete detaT;
 277   delete dphiT;
 278   delete cFlagT;
 279   delete cClusterT;
 280   delete idT;
 281   delete injet;
 282   delete sflag;
 283   delete hPtTotal;
 284   delete etCell;
 285   delete etaCell;
 286   delete phiCell;
 287   delete flagCell;
 288   delete etaJet;
 289   delete phiJet;
 290   delete etJet;
 291   delete etsigJet;
 292   delete etallJet;
 293   delete ncellsJet;
 294   delete multJet;
 295   delete idxjets;
 296   delete percentage;
 297   delete ncells;
 298   delete mult;
 299
 300
 301 }
 302
 303 ////////////////////////////////////////////////////////////////////////
 304
 305 void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell,
 306                                     Int_t* flagCell, Float_t etbgTotal, Double_t dEtTotal,
 307                                     Int_t& nJets, Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
 308                                     Float_t* etallJet, Int_t* ncellsJet)
 309 {
 310
 311   Int_t nCell = nIn;
 312
 313   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 314
 315   // Parameters from header
 316   Float_t minmove = header->GetMinMove();
 317   Float_t maxmove = header->GetMaxMove();
 318   Float_t rc      = header->GetRadius();
 319   Float_t etseed  = header->GetEtSeed();
 320
 321   // tmp array of jets form algoritm
 322   Float_t etaAlgoJet[30];
 323   Float_t phiAlgoJet[30];
 324   Float_t etAlgoJet[30];
 325   Int_t   ncellsAlgoJet[30];
 326
 327   //run algorithm//
 328
 329   // sort cells by et
 330   Int_t * index  = new Int_t[nCell];
 331   TMath::Sort(nCell, etCell, index);
 332
 333   // variable used in centroide loop
 334   Float_t eta   = 0.0;
 335   Float_t phi   = 0.0;
 336   Float_t eta0  = 0.0;
 337   Float_t phi0  = 0.0;
 338   Float_t etab  = 0.0;
 339   Float_t phib  = 0.0;
 340   Float_t etas  = 0.0;
 341   Float_t phis  = 0.0;
 342   Float_t ets   = 0.0;
 343   Float_t deta  = 0.0;
 344   Float_t dphi  = 0.0;
 345   Float_t dr    = 0.0;
 346   Float_t etsb  = 0.0;
 347   Float_t etasb = 0.0;
 348   Float_t phisb = 0.0;
 349   Float_t dphib = 0.0;
 350
 351
 352   for(Int_t icell = 0; icell < nCell; icell++){
 353       Int_t jcell = index[icell];
 354       if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
 355       if(flagCell[jcell] != 0) continue; // if cell was used before
 356       eta  = etaCell[jcell];
 357       phi  = phiCell[jcell];
 358       eta0 = eta;
 359       phi0 = phi;
 360       etab = eta;
 361       phib = phi;
 362       ets  = etCell[jcell];
 363       etas = 0.0;
 364       phis = 0.0;
 365       etsb = ets;
 366       etasb = 0.0;
 367       phisb = 0.0;
 368       for(Int_t kcell =0; kcell < nCell; kcell++){
 369           Int_t lcell = index[kcell];
 370           if(lcell == jcell) continue; // cell itself
 371           if(flagCell[lcell] != 0) continue; // cell used before
 372           if(etCell[lcell] > etCell[jcell]) continue;
 373           //calculate dr
 374           deta = etaCell[lcell] - eta;
 375           dphi = phiCell[lcell] - phi;
 376           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 377           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 378           dr = TMath::Sqrt(deta * deta + dphi * dphi);
 379           if(dr <= rc){
 380               // calculate offset from initiate cell
 381               deta = etaCell[lcell] - eta0;
 382               dphi = phiCell[lcell] - phi0;
 383               if (dphi < -TMath::Pi()) dphi = dphi + 2.0 * TMath::Pi();
 384               if (dphi >  TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
 385               etas = etas + etCell[lcell]*deta;
 386               phis = phis + etCell[lcell]*dphi;
 387               ets = ets + etCell[lcell];
 388                //new weighted eta and phi including this cell
 389                eta = eta0 + etas/ets;
 390                phi = phi0 + phis/ets;
 391                // if cone does not move much, just go to next step
 392                dphib = TMath::Abs(phi - phib);
 393                if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
 394                dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
 395                if(dr <= minmove) break;
 396                // cone should not move more than max_mov
 397                dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
 398                if(dr > maxmove){
 399                   eta = etab;
 400                   phi = phib;
 401                   ets = etsb;
 402                   etas = etasb;
 403                   phis = phisb;
 404                }else{ // store this loop information
 405                  etab=eta;
 406                  phib=phi;
 407                  etsb = ets;
 408                  etasb = etas;
 409                  phisb = phis;
 410                }
 411             }
 412         }//end of cells loop looking centroide
 413
 414         //avoid cones overloap (to be implemented in the future)
 415
 416         //flag cells in Rc, estimate total energy in cone
 417         Float_t etCone = 0.0;
 418         Int_t   nCellIn = 0;
 419         rc = header->GetRadius();
 420         for(Int_t ncell =0; ncell < nCell; ncell++){
 421             if(flagCell[ncell] != 0) continue; // cell used before
 422            //calculate dr
 423             deta = etaCell[ncell] - eta;
 424                  dphi = phiCell[ncell] - phi;
 425                  if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 426                  if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 427                  dr = TMath::Sqrt(deta * deta + dphi * dphi);
 428             if(dr <= rc){  // cell in cone
 429                flagCell[ncell] = -1;
 430                etCone+=etCell[ncell];
 431                nCellIn++;
 432             }
 433         }
 434
 435         // select jets with et > background
 436         // estimate max fluctuation of background in cone
 437         Double_t ncellin = (Double_t)nCellIn;
 438         Double_t ntcell  = (Double_t)nCell;
 439         Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
 440         // min cone et
 441         Double_t etcmin = etCone ;  // could be used etCone - etmin !!
 442         //desicions !! etbmax < etcmin
 443         for(Int_t mcell =0; mcell < nCell; mcell++){
 444             if(flagCell[mcell] == -1){
 445               if(etbmax < etcmin)
 446                  flagCell[mcell] = 1; //flag cell as used
 447               else
 448                  flagCell[mcell] = 0; // leave it free
 449             }
 450         }
 451         //store tmp jet info !!!
 452        if(etbmax < etcmin) {
 453              etaAlgoJet[nJets] = eta;
 454              phiAlgoJet[nJets] = phi;
 455              etAlgoJet[nJets] = etCone;
 456              ncellsAlgoJet[nJets] = nCellIn;
 457              nJets++;
 458         }
 459
 460   } // end of cells loop
 461
 462   //reorder jets by et in cone
 463   //sort jets by energy
 464   Int_t * idx  = new Int_t[nJets];
 465   TMath::Sort(nJets, etAlgoJet, idx);
 466   for(Int_t p = 0; p < nJets; p++){
 467      etaJet[p] = etaAlgoJet[idx[p]];
 468      phiJet[p] = phiAlgoJet[idx[p]];
 469      etJet[p] = etAlgoJet[idx[p]];
 470      etallJet[p] = etAlgoJet[idx[p]];
 471      ncellsJet[p] = ncellsAlgoJet[idx[p]];
 472   }
 473
 474
 475   //delete
 476   delete index;
 477   delete idx;
 478
 479 }
 480 ////////////////////////////////////////////////////////////////////////
 481
 482 void AliUA1JetFinderV2::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN,
 483                       Float_t* ptT, Float_t* etaT, Float_t* phiT,
 484                       Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
 485                       Int_t* multJet, Int_t* injet)
 486 {
 487   //background subtraction using cone method but without correction in dE/deta distribution
 488
 489   //calculate energy inside and outside cones
 490   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 491   Float_t rc= header->GetRadius();
 492   Float_t etIn[30];
 493   Float_t etOut = 0;
 494   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 495      // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
 496      for(Int_t ijet=0; ijet<nJ; ijet++){
 497          Float_t deta = etaT[jpart] - etaJet[ijet];
 498               Float_t dphi = phiT[jpart] - phiJet[ijet];
 499          if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 500               if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 501               Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 502          if(dr <= rc){ // particles inside this cone
 503              multJet[ijet]++;
 504              injet[jpart] = ijet;
 505              if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
 506                 etIn[ijet] += ptT[jpart];
 507                 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
 508              }
 509              break;
 510          }
 511      }// end jets loop
 512      if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
 513         etOut += ptT[jpart]; // particle outside cones and pt cut
 514   } //end particle loop
 515
 516   //estimate jets and background areas
 517   Float_t areaJet[30];
 518   Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
 519   for(Int_t k=0; k<nJ; k++){
 520       Float_t detamax = etaJet[k] + rc;
 521       Float_t detamin = etaJet[k] - rc;
 522       Float_t accmax = 0.0; Float_t accmin = 0.0;
 523       if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
 524          Float_t h = header->GetLegoEtaMax() - etaJet[k];
 525          accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
 526       }
 527       if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
 528          Float_t h = header->GetLegoEtaMax() + etaJet[k];
 529          accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
 530       }
 531       areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
 532       areaOut = areaOut - areaJet[k];
 533   }
 534   //subtract background using area method
 535   for(Int_t ljet=0; ljet<nJ; ljet++){
 536      Float_t areaRatio = areaJet[ljet]/areaOut;
 537      etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
 538   }
 539
 540   // estimate new total background
 541   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
 542   etbgTotalN = etOut*areaT/areaOut;
 543
 544
 545 }
 546
 547 ////////////////////////////////////////////////////////////////////////
 548
 549 void AliUA1JetFinderV2::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotalN,
 550                       Float_t* ptT, Float_t* etaT, Float_t* phiT,
 551                       Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
 552                       Int_t* multJet, Int_t* injet)
 553 {
 554
 555   //background subtraction using statistical method
 556   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 557   Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
 558
 559   //calculate energy inside
 560   Float_t rc= header->GetRadius();
 561   Float_t etIn[30];
 562
 563   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 564      //if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
 565      for(Int_t ijet=0; ijet<nJ; ijet++){
 566          Float_t deta = etaT[jpart] - etaJet[ijet];
 567               Float_t dphi = phiT[jpart] - phiJet[ijet];
 568          if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 569               if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 570               Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 571          if(dr <= rc){ // particles inside this cone
 572              multJet[ijet]++;
 573              injet[jpart] = ijet;
 574              if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
 575                 etIn[ijet]+= ptT[jpart];
 576                 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
 577              }
 578              break;
 579          }
 580      }// end jets loop
 581   } //end particle loop
 582
 583   //calc jets areas
 584   Float_t areaJet[30];
 585   Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
 586   for(Int_t k=0; k<nJ; k++){
 587       Float_t detamax = etaJet[k] + rc;
 588       Float_t detamin = etaJet[k] - rc;
 589       Float_t accmax = 0.0; Float_t accmin = 0.0;
 590       if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
 591          Float_t h = header->GetLegoEtaMax() - etaJet[k];
 592          accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
 593       }
 594       if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
 595          Float_t h = header->GetLegoEtaMax() + etaJet[k];
 596          accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
 597       }
 598       areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
 599   }
 600
 601   //subtract background using area method
 602   for(Int_t ljet=0; ljet<nJ; ljet++){
 603      Float_t areaRatio = areaJet[ljet]/areaOut;
 604      etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
 605   }
 606
 607   etbgTotalN = etbgStat;
 608
 609 }
 610
 611 ////////////////////////////////////////////////////////////////////////
 612
 613 void AliUA1JetFinderV2::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN,
 614                       Float_t* ptT, Float_t* etaT, Float_t* phiT,
 615                       Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
 616                       Int_t* multJet, Int_t* injet)
 617 {
 618    // Cone background subtraction method taking into acount dEt/deta distribution
 619     AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 620    //general
 621    Float_t rc= header->GetRadius();
 622    Float_t etamax = header->GetLegoEtaMax();
 623    Float_t etamin = header->GetLegoEtaMin();
 624    Int_t ndiv = 100;
 625
 626    // jet energy and area arrays
 627    TH1F* hEtJet[30];
 628    TH1F* hAreaJet[30];
 629    for(Int_t mjet=0; mjet<nJ; mjet++){
 630      char hEtname[256]; char hAreaname[256];
 631      sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
 632      hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
 633      hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
 634   }
 635    // background energy and area
 636    TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
 637    TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
 638
 639    //fill energies
 640    for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 641      for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
 642          Float_t deta = etaT[jpart] - etaJet[ijet];
 643          Float_t dphi = phiT[jpart] - phiJet[ijet];
 644          if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 645          if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 646          Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 647          if(dr <= rc){ // particles inside this cone
 648              injet[jpart] = ijet;
 649              multJet[ijet]++;
 650              if((fReader->GetCutFlag(jpart)) == 1){// pt cut
 651                 hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
 652                 if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
 653              }
 654              break;
 655          }
 656      }// end jets loop
 657      if(injet[jpart] == -1  && fReader->GetCutFlag(jpart) == 1)
 658         hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
 659   } //end particle loop
 660
 661    //calc areas
 662    Float_t eta0 = etamin;
 663    Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
 664    Float_t eta1 = eta0 + etaw;
 665    for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
 666       Float_t etac = eta0 + etaw/2.0;
 667       Float_t areabg = etaw*2.0*TMath::Pi();
 668       for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
 669           Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
 670           Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
 671           Float_t acc0 = 0.0; Float_t acc1 = 0.0;
 672           Float_t areaj = 0.0;
 673           if(deta0 > rc && deta1 < rc){
 674              acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
 675              areaj = acc1;
 676           }
 677           if(deta0 < rc && deta1 > rc){
 678              acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
 679              areaj = acc0;
 680           }
 681           if(deta0 < rc && deta1 < rc){
 682              acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
 683              acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
 684              if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
 685              if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
 686              if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
 687           }
 688           hAreaJet[ijet]->Fill(etac,areaj);
 689           areabg = areabg - areaj;
 690       } // end jets loop
 691       hAreaBackg->Fill(etac,areabg);
 692       eta0 = eta1;
 693       eta1 = eta1 + etaw;
 694    } // end loop for all eta bins
 695
 696    //subtract background
 697    for(Int_t kjet=0; kjet<nJ; kjet++){
 698        etJet[kjet] = 0.0; // first  clear etJet for this jet
 699        for(Int_t bin = 0; bin< ndiv; bin++){
 700            if(hAreaJet[kjet]->GetBinContent(bin)){
 701               Float_t areab = hAreaBackg->GetBinContent(bin);
 702               Float_t etb = hEtBackg->GetBinContent(bin);
 703               Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
 704               etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
 705            }
 706        }
 707    }
 708
 709    // calc background total
 710    Double_t etOut = hEtBackg->Integral();
 711    Double_t areaOut = hAreaBackg->Integral();
 712    Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
 713    etbgTotalN = etOut*areaT/areaOut;
 714
 715    //delete
 716    for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
 717        delete hEtJet[ljet];
 718        delete hAreaJet[ljet];
 719    }
 720
 721    delete hEtBackg;
 722    delete hAreaBackg;
 723 }
 724
 725 ////////////////////////////////////////////////////////////////////////
 726
 727
 728 void AliUA1JetFinderV2::SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN,
 729                       Float_t* ptT, Float_t* etaT, Float_t* phiT,
 730                       Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
 731                        Int_t* multJet, Int_t* injet)
 732 {
 733    // Ratio background subtraction method taking into acount dEt/deta distribution
 734     AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 735    //factor F calc before
 736     Float_t bgRatioCut = header->GetBackgCutRatio();
 737
 738
 739    //general
 740    Float_t rc= header->GetRadius();
 741    Float_t etamax = header->GetLegoEtaMax();
 742    Float_t etamin = header->GetLegoEtaMin();
 743    Int_t ndiv = 100;
 744
 745    // jet energy and area arrays
 746    TH1F* hEtJet[30];
 747    TH1F* hAreaJet[30];
 748    for(Int_t mjet=0; mjet<nJ; mjet++){
 749      char hEtname[256]; char hAreaname[256];
 750      sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
 751      hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);        // change range
 752      hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);  // change range
 753   }
 754    // background energy and area
 755    TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);         // change range
 756    TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);  // change range
 757
 758    //fill energies
 759    for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
 760      //if((fReader->GetCutFlag(jpart)) != 1) continue;
 761      for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
 762          Float_t deta = etaT[jpart] - etaJet[ijet];
 763               Float_t dphi = phiT[jpart] - phiJet[ijet];
 764          if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 765               if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 766               Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 767          if(dr <= rc){ // particles inside this cone
 768             multJet[ijet]++;
 769             injet[jpart] = ijet;
 770             if((fReader->GetCutFlag(jpart)) == 1){ //pt cut
 771                hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
 772                if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
 773             }
 774             break;
 775          }
 776      }// end jets loop
 777      if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
 778   } //end particle loop
 779
 780    //calc areas
 781    Float_t eta0 = etamin;
 782    Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
 783    Float_t eta1 = eta0 + etaw;
 784    for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
 785       Float_t etac = eta0 + etaw/2.0;
 786       Float_t areabg = etaw*2.0*TMath::Pi();
 787       for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
 788           Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
 789           Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
 790           Float_t acc0 = 0.0; Float_t acc1 = 0.0;
 791           Float_t areaj = 0.0;
 792           if(deta0 > rc && deta1 < rc){
 793              acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
 794              areaj = acc1;
 795           }
 796           if(deta0 < rc && deta1 > rc){
 797              acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
 798              areaj = acc0;
 799           }
 800           if(deta0 < rc && deta1 < rc){
 801              acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
 802              acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
 803              if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
 804              if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
 805              if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
 806           }
 807           hAreaJet[ijet]->Fill(etac,areaj);
 808           areabg = areabg - areaj;
 809       } // end jets loop
 810       hAreaBackg->Fill(etac,areabg);
 811       eta0 = eta1;
 812       eta1 = eta1 + etaw;
 813    } // end loop for all eta bins
 814
 815    //subtract background
 816    for(Int_t kjet=0; kjet<nJ; kjet++){
 817        etJet[kjet] = 0.0; // first  clear etJet for this jet
 818        for(Int_t bin = 0; bin< ndiv; bin++){
 819            if(hAreaJet[kjet]->GetBinContent(bin)){
 820               Float_t areab = hAreaBackg->GetBinContent(bin);
 821               Float_t etb = hEtBackg->GetBinContent(bin);
 822               Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
 823               etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
 824            }
 825        }
 826    }
 827
 828    // calc background total
 829    Double_t etOut = hEtBackg->Integral();
 830    Double_t areaOut = hAreaBackg->Integral();
 831    Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
 832    etbgTotalN = etOut*areaT/areaOut;
 833
 834    //delete
 835    for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
 836        delete hEtJet[ljet];
 837        delete hAreaJet[ljet];
 838    }
 839
 840    delete hEtBackg;
 841    delete hAreaBackg;
 842 }
 843
 844 ////////////////////////////////////////////////////////////////////////
 845
 846
 847 void AliUA1JetFinderV2::Reset()
 848 {
 849   fLego->Reset();
 850   fJets->ClearJets();
 851   AliJetFinder::Reset();
 852 }
 853
 854 ////////////////////////////////////////////////////////////////////////
 855
 856 void AliUA1JetFinderV2::WriteJHeaderToFile()
 857 {
 858   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 859   header->Write();
 860 }
 861
 862 ////////////////////////////////////////////////////////////////////////
 863
 864 void AliUA1JetFinderV2::Init()
 865 {
 866   // initializes some variables
 867   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 868    // book lego
 869   fLego = new
 870     TH2F("legoH","eta-phi",
 871          header->GetLegoNbinEta(), header->GetLegoEtaMin(),
 872          header->GetLegoEtaMax(),  header->GetLegoNbinPhi(),
 873          header->GetLegoPhiMin(),  header->GetLegoPhiMax());
 874
 875   fDebug = fReader->GetReaderHeader()->GetDebug();
 876   fOpt = fReader->GetReaderHeader()->GetDetector();
 877
 878   if(fOpt>0)
 879     fReader->CreateTasks();
 880
 881 }