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