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