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