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