JETAN/AliUA1JetFinderV2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //---------------------------------------------------------------------
  19 // UA1 Cone Algorithm Jet finder for charged + neutral jet studies
  20 // manages the search for jets using charged particle momentum and
  21 // neutral cell energy information
  22 // Based on UA1 V1 (from R. Diaz)
  23 // Author: magali.estienne@subatech.in2p3.fr
  24 //---------------------------------------------------------------------
  25
  26 #include <TClonesArray.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TLorentzVector.h>
  30 #include <TMath.h>
  31 #include <TRefArray.h>
  32 #include "TFile.h"
  33
  34 #include "AliUA1JetFinderV2.h"
  35 #include "AliUA1JetHeaderV1.h"
  36 #include "AliJetUnitArray.h"
  37
  38 class TArrayF;
  39 class TFile;
  40 class AliJetReader;
  41 class AliAODJet;
  42
  43 ClassImp(AliUA1JetFinderV2)
  44
  45
  46 ////////////////////////////////////////////////////////////////////////
  47 AliUA1JetFinderV2::AliUA1JetFinderV2() :
  48   AliJetFinder(),
  49   fLego(0),
  50   fDebug(0),
  51   fOpt(0)
  52 {
  53   //
  54   // Constructor
  55   //
  56 }
  57
  58 ////////////////////////////////////////////////////////////////////////
  59 AliUA1JetFinderV2::~AliUA1JetFinderV2()
  60 {
  61   //
  62   // Destructor
  63   //
  64 }
  65
  66 ////////////////////////////////////////////////////////////////////////
  67 void AliUA1JetFinderV2::FindJetsC()
  68 {
  69   //
  70   //  Used to find jets using charged particle momentum information
  71   //
  72   //  1) Fill cell map array
  73   //  2) calculate total energy and fluctuation level
  74   //  3) Run algorithm
  75   //     3.1) look centroides in cell map
  76   //     3.2) calculate total energy in cones
  77   //     3.3) flag as a possible jet
  78   //     3.4) reorder cones by energy
  79   //  4) subtract backg in accepted jets
  80   //  5) fill AliJet list
  81
  82   //  Transform input to pt,eta,phi plus lego
  83
  84   AliUA1JetHeaderV1* header  = (AliUA1JetHeaderV1*) fHeader;
  85   TClonesArray*      lvArray = fReader->GetMomentumArray();
  86   Int_t              nIn     = lvArray->GetEntries();
  87
  88   if (nIn == 0) return;
  89
  90   // local arrays for input
  91   Float_t* ptT    = new Float_t[nIn];
  92   Float_t* etaT   = new Float_t[nIn];
  93   Float_t* phiT   = new Float_t[nIn];
  94   Float_t* cFlagT = new Float_t[nIn]; // Temporarily added
  95   Float_t* sFlagT = new Float_t[nIn]; // Temporarily added
  96   Int_t*   injet  = new Int_t[nIn];
  97
  98   //total energy in array
  99   Float_t  etbgTotal = 0.0;
 100   TH1F*    hPtTotal  = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
 101
 102   // load input vectors and calculate total energy in array
 103   for (Int_t i = 0; i < nIn; i++){
 104     TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
 105     ptT[i]  = lv->Pt();
 106     etaT[i] = lv->Eta();
 107     phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
 108     cFlagT[i] = fReader->GetCutFlag(i);
 109     sFlagT[i] = fReader->GetSignalFlag(i);
 110
 111     if (fReader->GetCutFlag(i) != 1) continue;
 112     fLego->Fill(etaT[i], phiT[i], ptT[i]);
 113     hPtTotal->Fill(ptT[i]);
 114     etbgTotal+= ptT[i];
 115   }
 116
 117   fJets->SetNinput(nIn);
 118
 119   // calculate total energy and fluctuation in map
 120   Double_t meanpt   = hPtTotal->GetMean();
 121   Double_t ptRMS    = hPtTotal->GetRMS();
 122   Double_t npart    = hPtTotal->GetEntries();
 123   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
 124
 125   // arrays to hold jets
 126   Float_t* etaJet    = new Float_t[30];  // eta jet
 127   Float_t* phiJet    = new Float_t[30];  // phi jet
 128   Float_t* etJet     = new Float_t[30];  // et jet
 129   Float_t* etsigJet  = new Float_t[30];  // signal et in jet
 130   Float_t* etallJet  = new Float_t[30];  // total et in jet (tmp variable)
 131   Int_t*   ncellsJet = new Int_t[30];
 132   Int_t*   multJet   = new Int_t[30];
 133   //--- Added for jet reordering at the end of the jet finding procedure
 134   Float_t* etaJetOk    = new Float_t[30];
 135   Float_t* phiJetOk    = new Float_t[30];
 136   Float_t* etJetOk     = new Float_t[30];
 137   Float_t* etsigJetOk  = new Float_t[30];  // signal et in jet
 138   Float_t* etallJetOk  = new Float_t[30];  // total et in jet (tmp variable)
 139   Int_t*   ncellsJetOk = new Int_t[30];
 140   Int_t*   multJetOk   = new Int_t[30];
 141   //--------------------------
 142   Int_t nJets; // to hold number of jets found by algorithm
 143   Int_t nj;    // number of jets accepted
 144   Float_t prec  = header->GetPrecBg();
 145   Float_t bgprec = 1;
 146   while(bgprec > prec){
 147     //reset jet arrays in memory
 148     memset(etaJet,0,sizeof(Float_t)*30);
 149     memset(phiJet,0,sizeof(Float_t)*30);
 150     memset(etJet,0,sizeof(Float_t)*30);
 151     memset(etallJet,0,sizeof(Float_t)*30);
 152     memset(etsigJet,0,sizeof(Float_t)*30);
 153     memset(ncellsJet,0,sizeof(Int_t)*30);
 154     memset(multJet,0,sizeof(Int_t)*30);
 155     //--- Added for jet reordering at the end of the jet finding procedure
 156     memset(etaJetOk,0,sizeof(Float_t)*30);
 157     memset(phiJetOk,0,sizeof(Float_t)*30);
 158     memset(etJetOk,0,sizeof(Float_t)*30);
 159     memset(etallJetOk,0,sizeof(Float_t)*30);
 160     memset(etsigJetOk,0,sizeof(Float_t)*30);
 161     memset(ncellsJetOk,0,sizeof(Int_t)*30);
 162     memset(multJetOk,0,sizeof(Int_t)*30);
 163     //--------------------------
 164     nJets = 0;
 165     nj = 0;
 166
 167     // reset particles-jet array in memory
 168     memset(injet,-1,sizeof(Int_t)*nIn);
 169     //run cone algorithm finder
 170     RunAlgoritmC(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
 171
 172     //run background subtraction
 173     if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
 174       nj = header->GetNAcceptJets();
 175     else
 176       nj = nJets;
 177     //subtract background
 178     Float_t etbgTotalN = 0.0; //new background
 179     if(header->GetBackgMode() == 1) // standard
 180       SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 181     if(header->GetBackgMode() == 2) //cone
 182       SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 183     if(header->GetBackgMode() == 3) //ratio
 184       SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 185     if(header->GetBackgMode() == 4) //statistic
 186       SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 187     //calc precision
 188     if(etbgTotalN != 0.0)
 189       bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
 190     else
 191       bgprec = 0;
 192     etbgTotal = etbgTotalN; // update with new background estimation
 193   } //end while
 194
 195   // add jets to list
 196   Int_t* idxjets = new Int_t[nj];
 197   Int_t nselectj = 0;
 198   printf("Found %d jets \n", nj);
 199
 200   // Reorder jets by et in cone
 201   Int_t * idx  = new Int_t[nJets];
 202   TMath::Sort(nJets, etJet, idx);
 203   for(Int_t p = 0; p < nJets; p++){
 204     etaJetOk[p]    = etaJet[idx[p]];
 205     phiJetOk[p]    = phiJet[idx[p]];
 206     etJetOk[p]     = etJet[idx[p]];
 207     etallJetOk[p]  = etJet[idx[p]];
 208     etsigJetOk[p]  = etsigJet[idx[p]];
 209     ncellsJetOk[p] = ncellsJet[idx[p]];
 210     multJetOk[p]   = multJet[idx[p]];
 211   }
 212
 213   for(Int_t kj=0; kj<nj; kj++)
 214     {
 215       if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
 216           (etaJetOk[kj] < (header->GetJetEtaMin())) ||
 217           (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
 218       Float_t px, py,pz,en; // convert to 4-vector
 219       px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
 220       py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
 221       pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
 222       en = TMath::Sqrt(px * px + py * py + pz * pz);
 223       fJets->AddJet(px, py, pz, en);
 224
 225       AliAODJet jet(px, py, pz, en);
 226       jet.Print("");
 227
 228       AddJet(jet);
 229
 230       idxjets[nselectj] = kj;
 231       nselectj++;
 232     }
 233
 234   //add signal percentage and total signal  in AliJets for analysis tool
 235   Float_t* percentage  = new Float_t[nselectj];
 236   Int_t* ncells      = new Int_t[nselectj];
 237   Int_t* mult        = new Int_t[nselectj];
 238   for(Int_t i = 0; i< nselectj; i++)
 239     {
 240       percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
 241       ncells[i] = ncellsJetOk[idxjets[i]];
 242       mult[i] = multJetOk[idxjets[i]];
 243     }
 244
 245   //add particle-injet relationship ///
 246   for(Int_t bj = 0; bj < nIn; bj++)
 247     {
 248       if(injet[bj] == -1) continue; //background particle
 249       Int_t bflag = 0;
 250       for(Int_t ci = 0; ci< nselectj; ci++){
 251         if(injet[bj] == idxjets[ci]){
 252           injet[bj]= ci;
 253           bflag++;
 254           break;
 255         }
 256       }
 257       if(bflag == 0) injet[bj] = -1; // set as background particle
 258     }
 259
 260   fJets->SetNCells(ncells);
 261   fJets->SetPtFromSignal(percentage);
 262   fJets->SetMultiplicities(mult);
 263   fJets->SetInJet(injet);
 264   fJets->SetEtaIn(etaT);
 265   fJets->SetPhiIn(phiT);
 266   fJets->SetPtIn(ptT);
 267   fJets->SetEtAvg(etbgTotal/(4*(header->GetLegoEtaMax())*TMath::Pi()));
 268
 269   //delete
 270   delete[] ptT;
 271   delete[] etaT;
 272   delete[] phiT;
 273   delete[] cFlagT;
 274   delete[] sFlagT;
 275   delete[] injet;
 276   delete[] hPtTotal;
 277   delete[] etaJet;
 278   delete[] phiJet;
 279   delete[] etJet;
 280   delete[] etsigJet;
 281   delete[] etallJet;
 282   delete[] ncellsJet;
 283   delete[] multJet;
 284   delete[] idxjets;
 285   delete[] percentage;
 286   delete[] ncells;
 287   delete[] mult;
 288   //--- Added for jet reordering
 289   delete etaJetOk;
 290   delete phiJetOk;
 291   delete etJetOk;
 292   delete etsigJetOk;
 293   delete etallJetOk;
 294   delete ncellsJetOk;
 295   delete multJetOk;
 296   //--------------------------
 297
 298 }
 299
 300 ////////////////////////////////////////////////////////////////////////
 301 void AliUA1JetFinderV2::FindJets()
 302 {
 303   //
 304   //  Used to find jets using charged particle momentum information
 305   //  & neutral energy from calo cells
 306   //
 307   //  1) Fill cell map array
 308   //  2) calculate total energy and fluctuation level
 309   //  3) Run algorithm
 310   //     3.1) look centroides in cell map
 311   //     3.2) calculate total energy in cones
 312   //     3.3) flag as a possible jet
 313   //     3.4) reorder cones by energy
 314   //  4) subtract backg in accepted jets
 315   //  5) fill AliJet list
 316
 317   // transform input to pt,eta,phi plus lego
 318
 319   AliUA1JetHeaderV1* header   = (AliUA1JetHeaderV1*) fHeader;
 320   TClonesArray*      fUnit    = fReader->GetUnitArray();
 321   Int_t              nCand    = fReader->GetNumCandidate();
 322   Int_t              nCandCut = fReader->GetNumCandidateCut();
 323   Int_t              nIn      = fUnit->GetEntries();
 324   Float_t            ptMin   = fReader->GetReaderHeader()->GetPtCut();
 325
 326   fDebug   = fReader->GetReaderHeader()->GetDebug();
 327   if (nIn == 0) return;
 328
 329   Int_t nCandidateCut = 0;
 330   Int_t nCandidate = 0;
 331
 332   nCandidate = nCand;
 333   nCandidateCut = nCandCut;
 334
 335   // local arrays for input No Cuts
 336   // Both pt < ptMin and pt > ptMin
 337   Float_t*   ptT       = new Float_t[nCandidate];
 338   Float_t*   en2T      = new Float_t[nCandidate];
 339   Float_t*   pt2T      = new Float_t[nCandidate];
 340   Int_t*     detT      = new Int_t[nCandidate];
 341   Float_t*   etaT      = new Float_t[nCandidate];
 342   Float_t*   phiT      = new Float_t[nCandidate];
 343   Float_t*   cFlagT    = new Float_t[nCandidate];
 344   Float_t*   cFlag2T   = new Float_t[nCandidate];
 345   Float_t*   sFlagT    = new Float_t[nCandidate];
 346   Float_t*   cClusterT = new Float_t[nCandidate];
 347   Int_t*     vectT     = new Int_t[nCandidate];
 348   Int_t      loop1     = 0;
 349   Int_t*     injet     = new Int_t[nCandidate];
 350   Int_t*     sflag     = new Int_t[nCandidate];
 351   TRefArray* trackRef  = new TRefArray();
 352
 353   //total energy in array
 354   Float_t  etbgTotal = 0.0;
 355   TH1F*    hPtTotal  = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
 356
 357   // Input cell info
 358   Float_t *etCell    = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
 359   Float_t *etaCell   = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
 360   Float_t *phiCell   = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
 361   Int_t   *flagCell  = new Int_t[nIn];   //! Cell phi - Extracted from UnitArray
 362   Float_t *etCell2   = new Float_t[nIn];  //! Cell Energy - Extracted from UnitArray
 363   Float_t *etaCell2  = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
 364   Float_t *phiCell2  = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
 365   Int_t   *flagCell2 = new Int_t[nIn];  //! Cell phi - Extracted from UnitArray
 366
 367   // Information extracted from fUnitArray
 368   // Load input vectors and calculate total energy in array
 369   for(Int_t i=0; i<nIn; i++)
 370     {
 371       // Recover particle information from UnitArray
 372
 373       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
 374       TRefArray* ref = uArray->GetUnitTrackRef();
 375       Int_t nRef = ref->GetEntries();
 376
 377       if(uArray->GetUnitEnergy()>0.){
 378
 379         for(Int_t jpart=0; jpart<nRef;jpart++)
 380           trackRef->Add((AliVTrack*)ref->At(jpart));
 381         ptT[loop1]   = uArray->GetUnitEnergy();
 382         detT[loop1]  = uArray->GetUnitDetectorFlag();
 383         etaT[loop1]  = uArray->GetUnitEta();
 384         phiT[loop1]  = uArray->GetUnitPhi();
 385         cFlagT[loop1]= uArray->GetUnitCutFlag();   // pt cut tpc
 386         cFlag2T[loop1]= uArray->GetUnitCutFlag2(); // pt cut emcal
 387         sFlagT[loop1]= uArray->GetUnitSignalFlag();
 388         vectT[loop1] = nRef;
 389         if(cFlagT[loop1] == 1 || cFlag2T[loop1] == 1) {
 390           pt2T[loop1] = 0.;
 391           en2T[loop1] = 0.;
 392           if(detT[loop1]==1){
 393             en2T[loop1] = ptT[loop1] - header->GetMinCellEt();
 394             if(en2T[loop1] < 0) en2T[loop1]=0;
 395             hPtTotal->Fill(en2T[loop1]);
 396             etbgTotal += en2T[loop1];
 397           }
 398           if(detT[loop1]==0){ // TPC+ITS
 399             Float_t pt = 0.;
 400             for(Int_t j=0; j<nRef;j++){
 401               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
 402               x = ((AliVTrack*)ref->At(j))->Px();
 403               y = ((AliVTrack*)ref->At(j))->Py();
 404               z = ((AliVTrack*)ref->At(j))->Pz();
 405               pt = TMath::Sqrt(x*x+y*y);
 406               if(pt>ptMin) {
 407                 pt2T[loop1] += pt;
 408                 en2T[loop1] += pt;
 409                 hPtTotal->Fill(pt);
 410                 etbgTotal+= pt;
 411               }
 412             }
 413           }
 414           if(detT[loop1]==2) { // EMCal
 415             Float_t ptCTot = 0.;
 416             Float_t pt = 0.;
 417             Float_t enC = 0.;
 418             for(Int_t j=0; j<nRef;j++){
 419               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
 420               x = ((AliVTrack*)ref->At(j))->Px();
 421               y = ((AliVTrack*)ref->At(j))->Py();
 422               z = ((AliVTrack*)ref->At(j))->Pz();
 423               pt = TMath::Sqrt(x*x+y*y);
 424               if(pt>ptMin) {
 425                 pt2T[loop1]+=pt;
 426                 en2T[loop1]+=pt;
 427                 hPtTotal->Fill(pt);
 428                 etbgTotal+= pt;
 429               }
 430               ptCTot += pt;
 431             }
 432             enC = ptT[loop1] - ptCTot - header->GetMinCellEt();
 433             if(enC < 0.) enC=0.;
 434             en2T[loop1] += enC;
 435             hPtTotal->Fill(enC);
 436             etbgTotal+= enC;
 437           }
 438         }
 439         loop1++;
 440       }
 441
 442       if(uArray->GetUnitCutFlag()==1) {
 443         if(uArray->GetUnitDetectorFlag()==1){ // EMCal case
 444           etCell[i] = uArray->GetUnitEnergy() - header->GetMinCellEt();
 445           if ((uArray->GetUnitEnergy() - header->GetMinCellEt()) < 0.0) etCell[i]=0.;
 446           etaCell[i] = uArray->GetUnitEta();
 447           phiCell[i] = uArray->GetUnitPhi();
 448           flagCell[i] = 0; // default
 449           etCell2[i] = etCell[i];
 450           etaCell2[i] = uArray->GetUnitEta();
 451           phiCell2[i] = uArray->GetUnitPhi();
 452           flagCell2[i] = 0; // default
 453         }
 454         if(uArray->GetUnitDetectorFlag()==0){ // TPC case
 455           Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
 456           for(Int_t j=0; j<nRef;j++)
 457             {
 458               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
 459               x = ((AliVTrack*)ref->At(j))->Px();
 460               y = ((AliVTrack*)ref->At(j))->Py();
 461               z = ((AliVTrack*)ref->At(j))->Pz();
 462               pt = TMath::Sqrt(x*x+y*y);
 463               if(pt>ptMin) {
 464                 et1 += pt;
 465                 et2 += pt;
 466               }
 467             }
 468           etCell[i] = et1;
 469           etCell2[i] = et2;
 470           if(et1 < 0.) etCell[i] = etCell2[i] = 0.;
 471           etaCell[i] = uArray->GetUnitEta();
 472           phiCell[i] = uArray->GetUnitPhi();
 473           flagCell[i] = 0; // default
 474           etaCell2[i] = uArray->GetUnitEta();
 475           phiCell2[i] = uArray->GetUnitPhi();
 476           flagCell2[i] = 0; // default
 477         }
 478         if(uArray->GetUnitDetectorFlag()==2){ // TPC + EMCal case
 479           Float_t ptCTot = 0.;
 480           Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
 481           Float_t enC = 0.;
 482           for(Int_t j=0; j<nRef;j++)
 483             {
 484               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
 485               x = ((AliVTrack*)ref->At(j))->Px();
 486               y = ((AliVTrack*)ref->At(j))->Py();
 487               z = ((AliVTrack*)ref->At(j))->Pz();
 488               pt = TMath::Sqrt(x*x+y*y);
 489               if(pt>ptMin) {
 490                 et1 += pt;
 491                 et2 += pt;
 492               }
 493               ptCTot += pt;
 494             }
 495           enC = uArray->GetUnitEnergy() - ptCTot;
 496           etCell[i] = et1 + enC - header->GetMinCellEt();
 497           etCell2[i] = et2 + enC - header->GetMinCellEt();
 498           if((enC + et1 - header->GetMinCellEt()) < 0.) etCell[i] = etCell2[i] = 0.;
 499           etaCell[i] = uArray->GetUnitEta();
 500           phiCell[i] = uArray->GetUnitPhi();
 501           flagCell[i] = 0; // default
 502           etaCell2[i] = uArray->GetUnitEta();
 503           phiCell2[i] = uArray->GetUnitPhi();
 504           flagCell2[i] = 0; // default
 505         }
 506       }
 507       else {
 508         etCell[i]  = 0.;
 509         etaCell[i] = uArray->GetUnitEta();
 510         phiCell[i] = uArray->GetUnitPhi();
 511         flagCell[i] = 0;
 512         etCell2[i]  = 0.;
 513         etaCell2[i] = uArray->GetUnitEta();
 514         phiCell2[i] = uArray->GetUnitPhi();
 515         flagCell2[i] = 0;
 516       }
 517     } // end loop on nCandidate
 518
 519   fJets->SetNinput(nCandidate);
 520
 521   // calculate total energy and fluctuation in map
 522   Double_t meanpt = hPtTotal->GetMean();
 523   Double_t ptRMS = hPtTotal->GetRMS();
 524   Double_t npart = hPtTotal->GetEntries();
 525   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
 526
 527   // arrays to hold jets
 528   Float_t* etaJet    = new Float_t[30];
 529   Float_t* phiJet    = new Float_t[30];
 530   Float_t* etJet     = new Float_t[30];
 531   Float_t* etsigJet  = new Float_t[30]; //signal et in jet
 532   Float_t* etallJet  = new Float_t[30];  // total et in jet (tmp variable)
 533   Int_t*   ncellsJet = new Int_t[30];
 534   Int_t*   multJet   = new Int_t[30];
 535   //--- Added by me for jet reordering at the end of the jet finding procedure
 536   Float_t* etaJetOk    = new Float_t[30];
 537   Float_t* phiJetOk    = new Float_t[30];
 538   Float_t* etJetOk     = new Float_t[30];
 539   Float_t* etsigJetOk  = new Float_t[30]; //signal et in jet
 540   Float_t* etallJetOk  = new Float_t[30];  // total et in jet (tmp variable)
 541   Int_t*   ncellsJetOk = new Int_t[30];
 542   Int_t*   multJetOk   = new Int_t[30];
 543   //--------------------------
 544   Int_t    nJets; // to hold number of jets found by algorithm
 545   Int_t    nj;    // number of jets accepted
 546   Float_t  prec  = header->GetPrecBg();
 547   Float_t  bgprec = 1;
 548
 549   while(bgprec > prec){
 550
 551     //reset jet arrays in memory
 552     memset(etaJet,0,sizeof(Float_t)*30);
 553     memset(phiJet,0,sizeof(Float_t)*30);
 554     memset(etJet,0,sizeof(Float_t)*30);
 555     memset(etallJet,0,sizeof(Float_t)*30);
 556     memset(etsigJet,0,sizeof(Float_t)*30);
 557     memset(ncellsJet,0,sizeof(Int_t)*30);
 558     memset(multJet,0,sizeof(Int_t)*30);
 559     //--- Added by me for jet reordering at the end of the jet finding procedure
 560     memset(etaJetOk,0,sizeof(Float_t)*30);
 561     memset(phiJetOk,0,sizeof(Float_t)*30);
 562     memset(etJetOk,0,sizeof(Float_t)*30);
 563     memset(etallJetOk,0,sizeof(Float_t)*30);
 564     memset(etsigJetOk,0,sizeof(Float_t)*30);
 565     memset(ncellsJetOk,0,sizeof(Int_t)*30);
 566     memset(multJetOk,0,sizeof(Int_t)*30);
 567
 568     nJets = 0;
 569     nj = 0;
 570
 571     // reset particles-jet array in memory
 572     memset(injet,-1,sizeof(Int_t)*nCandidate);
 573     //run cone algorithm finder
 574     RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etCell2,etaCell2,phiCell2,
 575                 flagCell2,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,
 576                 etallJet,ncellsJet);
 577
 578     //run background subtraction
 579     if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
 580       nj = header->GetNAcceptJets();
 581     else
 582       nj = nJets;
 583
 584     //subtract background
 585     Float_t etbgTotalN = 0.0; //new background
 586     if(header->GetBackgMode() == 1) // standard
 587       SubtractBackg(nCandidate,nj,etbgTotalN,en2T,vectT,etaT,phiT,cFlagT,cFlag2T,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 588     // To be modified ------------------------
 589     if(header->GetBackgMode() == 2) //cone
 590       SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 591     if(header->GetBackgMode() == 3) //ratio
 592       SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 593     if(header->GetBackgMode() == 4) //statistic
 594       SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 595     //----------------------------------------
 596     //calc precision
 597     if(etbgTotalN != 0.0)
 598       bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
 599     else
 600       bgprec = 0;
 601     etbgTotal = etbgTotalN; // update with new background estimation
 602   } //end while
 603
 604   // add jets to list
 605   Int_t* idxjets = new Int_t[nj];
 606   Int_t nselectj = 0;
 607   printf("Found %d jets \n", nj);
 608
 609   // Reorder jets by et in cone
 610   // Sort jets by energy
 611   Int_t * idx  = new Int_t[nJets];
 612   TMath::Sort(nJets, etJet, idx);
 613   for(Int_t p = 0; p < nJets; p++)
 614     {
 615       etaJetOk[p]    = etaJet[idx[p]];
 616       phiJetOk[p]    = phiJet[idx[p]];
 617       etJetOk[p]     = etJet[idx[p]];
 618       etallJetOk[p]  = etJet[idx[p]];
 619       etsigJetOk[p]  = etsigJet[idx[p]];
 620       ncellsJetOk[p] = ncellsJet[idx[p]];
 621       multJetOk[p]   = multJet[idx[p]];
 622     }
 623
 624   for(Int_t kj=0; kj<nj; kj++)
 625     {
 626       if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
 627           (etaJetOk[kj] < (header->GetJetEtaMin())) ||
 628           (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
 629       Float_t px, py,pz,en; // convert to 4-vector
 630       px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
 631       py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
 632       pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
 633       en = TMath::Sqrt(px * px + py * py + pz * pz);
 634       fJets->AddJet(px, py, pz, en);
 635       AliAODJet jet(px, py, pz, en);
 636       jet.Print("");
 637
 638       AddJet(jet);
 639
 640       idxjets[nselectj] = kj;
 641       nselectj++;
 642     }
 643
 644   //add signal percentage and total signal  in AliJets for analysis tool
 645   Float_t* percentage  = new Float_t[nselectj];
 646   Int_t* ncells      = new Int_t[nselectj];
 647   Int_t* mult        = new Int_t[nselectj];
 648   for(Int_t i = 0; i< nselectj; i++)
 649     {
 650       percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
 651       ncells[i]     = ncellsJetOk[idxjets[i]];
 652       mult[i]       = multJetOk[idxjets[i]];
 653     }
 654
 655   //add particle-injet relationship ///
 656   for(Int_t bj = 0; bj < nCandidate; bj++)
 657     {
 658       if(injet[bj] == -1) continue; //background particle
 659       Int_t bflag = 0;
 660       for(Int_t ci = 0; ci< nselectj; ci++){
 661         if(injet[bj] == idxjets[ci]){
 662           injet[bj]= ci;
 663           bflag++;
 664           break;
 665         }
 666       }
 667     if(bflag == 0) injet[bj] = -1; // set as background particle
 668   }
 669
 670   fJets->SetNCells(ncells);
 671   fJets->SetPtFromSignal(percentage);
 672   fJets->SetMultiplicities(mult);
 673   fJets->SetInJet(injet);
 674   fJets->SetEtaIn(etaT);
 675   fJets->SetPhiIn(phiT);
 676   fJets->SetPtIn(ptT);
 677   fJets->SetTrackReferences(trackRef);
 678   fJets->SetDetectorFlagIn(detT);
 679   fJets->SetEtAvg(etbgTotal/(2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax()-header->GetLegoPhiMin())));
 680
 681   //delete
 682   delete ptT;
 683   delete en2T;
 684   delete pt2T;
 685   delete etaT;
 686   delete phiT;
 687   trackRef->Delete();
 688   delete trackRef;
 689   delete detT;
 690   delete cFlagT;
 691   delete cFlag2T;
 692   delete sFlagT;
 693   delete cClusterT;
 694   delete vectT;
 695   delete injet;
 696   delete sflag;
 697   delete hPtTotal;
 698   delete etCell;
 699   delete etaCell;
 700   delete phiCell;
 701   delete flagCell;
 702   delete etCell2;
 703   delete etaCell2;
 704   delete phiCell2;
 705   delete flagCell2;
 706   delete etaJet;
 707   delete phiJet;
 708   delete etJet;
 709   delete etsigJet;
 710   delete etallJet;
 711   delete ncellsJet;
 712   delete multJet;
 713   //--- Added for jet reordering
 714   delete etaJetOk;
 715   delete phiJetOk;
 716   delete etJetOk;
 717   delete etsigJetOk;
 718   delete etallJetOk;
 719   delete ncellsJetOk;
 720   delete multJetOk;
 721   //--------------------------
 722   delete idxjets;
 723   delete percentage;
 724   delete ncells;
 725   delete mult;
 726
 727 }
 728
 729 ////////////////////////////////////////////////////////////////////////
 730 void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell,
 731                                     Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2,
 732                                     const Float_t* phiCell2, const Int_t* flagCell2, Float_t etbgTotal,
 733                                     Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet,
 734                                     Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet)
 735 {
 736   //
 737   // Main method for jet finding
 738   // UA1 base cone finder
 739   //
 740
 741   Int_t nCell  = nIn;
 742   fDebug = fReader->GetReaderHeader()->GetDebug();
 743
 744   // Dump lego
 745   // Check enough space! *to be done*
 746   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 747   for(Int_t i=0; i<nCell; i++){
 748     etCell[i]   = etCell2[i];
 749     etaCell[i]  = etaCell2[i];
 750     phiCell[i]  = phiCell2[i];
 751     flagCell[i] = flagCell2[i];
 752   }
 753
 754   // Parameters from header
 755   Float_t minmove = header->GetMinMove();
 756   Float_t maxmove = header->GetMaxMove();
 757   Float_t rc      = header->GetRadius();
 758   Float_t etseed  = header->GetEtSeed();
 759
 760   // Tmp array of jets form algoritm
 761   Float_t etaAlgoJet[30];
 762   Float_t phiAlgoJet[30];
 763   Float_t etAlgoJet[30];
 764   Int_t   ncellsAlgoJet[30];
 765
 766   // Run algorithm//
 767
 768   // Sort cells by et
 769   Int_t * index  = new Int_t[nCell];
 770   TMath::Sort(nCell, etCell, index);
 771
 772   // Variable used in centroide loop
 773   Float_t eta   = 0.0;
 774   Float_t phi   = 0.0;
 775   Float_t eta0  = 0.0;
 776   Float_t phi0  = 0.0;
 777   Float_t etab  = 0.0;
 778   Float_t phib  = 0.0;
 779   Float_t etas  = 0.0;
 780   Float_t phis  = 0.0;
 781   Float_t ets   = 0.0;
 782   Float_t deta  = 0.0;
 783   Float_t dphi  = 0.0;
 784   Float_t dr    = 0.0;
 785   Float_t etsb  = 0.0;
 786   Float_t etasb = 0.0;
 787   Float_t phisb = 0.0;
 788   Float_t dphib = 0.0;
 789
 790   for(Int_t icell = 0; icell < nCell; icell++)
 791     {
 792       Int_t jcell = index[icell];
 793       if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
 794       if(flagCell[jcell] != 0) continue; // if cell was used before
 795
 796       eta  = etaCell[jcell];
 797       phi  = phiCell[jcell];
 798       eta0 = eta;
 799       phi0 = phi;
 800       etab = eta;
 801       phib = phi;
 802       ets  = etCell[jcell];
 803       etas = 0.0;
 804       phis = 0.0;
 805       etsb = ets;
 806       etasb = 0.0;
 807       phisb = 0.0;
 808       for(Int_t kcell =0; kcell < nCell; kcell++)
 809         {
 810           Int_t lcell = index[kcell];
 811           if(lcell == jcell) continue; // cell itself
 812           if(flagCell[lcell] != 0) continue; // cell used before
 813           if(etCell[lcell] > etCell[jcell]) continue;  // can this happen
 814           //calculate dr
 815           deta = etaCell[lcell] - eta;
 816           dphi = TMath::Abs(phiCell[lcell] - phi);
 817           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 818           dr = TMath::Sqrt(deta * deta + dphi * dphi);
 819           if(dr <= rc){
 820             // calculate offset from initiate cell
 821             deta = etaCell[lcell] - eta0;
 822             dphi = phiCell[lcell] - phi0;
 823             if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
 824             if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
 825             etas = etas + etCell[lcell]*deta;
 826             phis = phis + etCell[lcell]*dphi;
 827             ets = ets + etCell[lcell];
 828             //new weighted eta and phi including this cell
 829             eta = eta0 + etas/ets;
 830             phi = phi0 + phis/ets;
 831             // if cone does not move much, just go to next step
 832             dphib = TMath::Abs(phi - phib);
 833             if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
 834             dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
 835             if(dr <= minmove) break;
 836             // cone should not move more than max_mov
 837             dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
 838             if(dr > maxmove){
 839               eta = etab;
 840               phi = phib;
 841               ets = etsb;
 842               etas = etasb;
 843               phis = phisb;
 844             } else { // store this loop information
 845               etab  = eta;
 846               phib  = phi;
 847               etsb  = ets;
 848               etasb = etas;
 849               phisb = phis;
 850             }
 851           } // inside cone
 852         }//end of cells loop looking centroide
 853
 854         //avoid cones overloap (to be implemented in the future)
 855
 856         //flag cells in Rc, estimate total energy in cone
 857       Float_t etCone = 0.0;
 858       Int_t   nCellIn = 0;
 859       Int_t   nCellOut = 0;
 860       rc = header->GetRadius();
 861
 862       for(Int_t ncell =0; ncell < nCell; ncell++)
 863         {
 864           if(flagCell[ncell] != 0) continue; // cell used before
 865           //calculate dr
 866           deta = etaCell[ncell] - eta;
 867           //        if(deta <= rc){ // Added to improve velocity -> to be tested
 868           dphi = phiCell[ncell] - phi;
 869           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 870           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 871           //          if(dphi <= rc){ // Added to improve velocity -> to be tested
 872           dr = TMath::Sqrt(deta * deta + dphi * dphi);
 873           if(dr <= rc){  // cell in cone
 874             flagCell[ncell] = -1;
 875             etCone+=etCell[ncell];
 876             nCellIn++;
 877           }
 878           else nCellOut++;
 879           //    } // end deta <= rc
 880           //        } // end dphi <= rc
 881         }
 882
 883       // select jets with et > background
 884       // estimate max fluctuation of background in cone
 885       Double_t ncellin = (Double_t)nCellIn;
 886       Double_t ntcell  = (Double_t)nCell;
 887       Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/(ntcell));
 888       // min cone et
 889       Double_t etcmin = etCone ;  // could be used etCone - etmin !!
 890       //decisions !! etbmax < etcmin
 891
 892       for(Int_t mcell =0; mcell < nCell; mcell++)
 893         {
 894           if(flagCell[mcell] == -1){
 895             if(etbmax < etcmin)
 896               flagCell[mcell] = 1; //flag cell as used
 897             else
 898               flagCell[mcell] = 0; // leave it free
 899           }
 900         }
 901       //store tmp jet info !!!
 902       if(etbmax < etcmin)
 903         {
 904           etaAlgoJet[nJets]    = eta;
 905           phiAlgoJet[nJets]    = phi;
 906           etAlgoJet[nJets]     = etCone;
 907           ncellsAlgoJet[nJets] = nCellIn;
 908           nJets++;
 909         }
 910
 911     } // end of cells loop
 912
 913   for(Int_t p = 0; p < nJets; p++)
 914     {
 915       etaJet[p]    = etaAlgoJet[p];
 916       phiJet[p]    = phiAlgoJet[p];
 917       etJet[p]     = etAlgoJet[p];
 918       etallJet[p]  = etAlgoJet[p];
 919       ncellsJet[p] = ncellsAlgoJet[p];
 920     }
 921
 922   //delete
 923   delete index;
 924
 925 }
 926
 927 ////////////////////////////////////////////////////////////////////////
 928 void AliUA1JetFinderV2::RunAlgoritmC(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
 929                                   Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
 930                                   Float_t* etallJet, Int_t* ncellsJet)
 931 {
 932   // Dump lego
 933   // Check enough space! *to be done*
 934   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
 935   Float_t etCell[60000];   //! Cell Energy
 936   Float_t etaCell[60000];  //! Cell eta
 937   Float_t phiCell[60000];  //! Cell phi
 938   Int_t   flagCell[60000]; //! Cell flag
 939
 940   Int_t nCell = 0;
 941   TAxis* xaxis = fLego->GetXaxis();
 942   TAxis* yaxis = fLego->GetYaxis();
 943   Float_t e = 0.0;
 944   for (Int_t i = 1; i <= header->GetLegoNbinEta(); i++)
 945     {
 946       for (Int_t j = 1; j <= header->GetLegoNbinPhi(); j++)
 947         {
 948           e = fLego->GetBinContent(i,j);
 949           if (e < 0.0) continue; // don't include this cells
 950           Float_t eta  = xaxis->GetBinCenter(i);
 951           Float_t phi  = yaxis->GetBinCenter(j);
 952           etCell[nCell]  = e;
 953           etaCell[nCell] = eta;
 954           phiCell[nCell] = phi;
 955           flagCell[nCell] = 0; //default
 956           nCell++;
 957         }
 958     }
 959
 960   // Parameters from header
 961   Float_t minmove = header->GetMinMove();
 962   Float_t maxmove = header->GetMaxMove();
 963   Float_t rc      = header->GetRadius();
 964   Float_t etseed  = header->GetEtSeed();
 965
 966   // Tmp array of jets form algoritm
 967   Float_t etaAlgoJet[30];
 968   Float_t phiAlgoJet[30];
 969   Float_t etAlgoJet[30];
 970   Int_t   ncellsAlgoJet[30];
 971
 972   // Run algorithm//
 973
 974   // Sort cells by et
 975   Int_t * index  = new Int_t[nCell];
 976   TMath::Sort(nCell, etCell, index);
 977   // variable used in centroide loop
 978   Float_t eta   = 0.0;
 979   Float_t phi   = 0.0;
 980   Float_t eta0  = 0.0;
 981   Float_t phi0  = 0.0;
 982   Float_t etab  = 0.0;
 983   Float_t phib  = 0.0;
 984   Float_t etas  = 0.0;
 985   Float_t phis  = 0.0;
 986   Float_t ets   = 0.0;
 987   Float_t deta  = 0.0;
 988   Float_t dphi  = 0.0;
 989   Float_t dr    = 0.0;
 990   Float_t etsb  = 0.0;
 991   Float_t etasb = 0.0;
 992   Float_t phisb = 0.0;
 993   Float_t dphib = 0.0;
 994
 995   for(Int_t icell = 0; icell < nCell; icell++)
 996     {
 997       Int_t jcell = index[icell];
 998       if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
 999       if(flagCell[jcell] != 0) continue; // if cell was used before
1000
1001       eta  = etaCell[jcell];
1002       phi  = phiCell[jcell];
1003       eta0 = eta;
1004       phi0 = phi;
1005       etab = eta;
1006       phib = phi;
1007       ets  = etCell[jcell];
1008       etas = 0.0;
1009       phis = 0.0;
1010       etsb = ets;
1011       etasb = 0.0;
1012       phisb = 0.0;
1013       for(Int_t kcell =0; kcell < nCell; kcell++)
1014         {
1015           Int_t lcell = index[kcell];
1016           if(lcell == jcell) continue; // cell itself
1017           if(flagCell[lcell] != 0) continue; // cell used before
1018           if(etCell[lcell] > etCell[jcell]) continue; // can this happen
1019           //calculate dr
1020           deta = etaCell[lcell] - eta;
1021           dphi = TMath::Abs(phiCell[lcell] - phi);
1022           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1023           dr = TMath::Sqrt(deta * deta + dphi * dphi);
1024           if(dr <= rc)
1025             {
1026               // calculate offset from initiate cell
1027               deta = etaCell[lcell] - eta0;
1028               dphi = phiCell[lcell] - phi0;
1029               if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
1030               if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
1031               etas = etas + etCell[lcell]*deta;
1032               phis = phis + etCell[lcell]*dphi;
1033               ets = ets + etCell[lcell];
1034               //new weighted eta and phi including this cell
1035               eta = eta0 + etas/ets;
1036               phi = phi0 + phis/ets;
1037               // if cone does not move much, just go to next step
1038               dphib = TMath::Abs(phi - phib);
1039               if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
1040               dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
1041               if(dr <= minmove) break;
1042               // cone should not move more than max_mov
1043               dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
1044               if(dr > maxmove){
1045                 eta = etab;
1046                 phi = phib;
1047                 ets = etsb;
1048                 etas = etasb;
1049                 phis = phisb;
1050               } else { // store this loop information
1051                 etab=eta;
1052                 phib=phi;
1053                 etsb = ets;
1054                 etasb = etas;
1055                 phisb = phis;
1056                 }
1057             } // inside cone
1058         }//end of cells loop looking centroide
1059
1060       // Avoid cones overloap (to be implemented in the future)
1061
1062       // Flag cells in Rc, estimate total energy in cone
1063       Float_t etCone   = 0.0;
1064       Int_t   nCellIn  = 0;
1065       Int_t   nCellOut = 0;
1066       rc = header->GetRadius();
1067       for(Int_t ncell =0; ncell < nCell; ncell++)
1068         {
1069           if(flagCell[ncell] != 0) continue; // cell used before
1070           //calculate dr
1071           deta = etaCell[ncell] - eta;
1072           dphi = phiCell[ncell] - phi;
1073           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1074           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1075           dr = TMath::Sqrt(deta * deta + dphi * dphi);
1076           if(dr <= rc){  // cell in cone
1077             flagCell[ncell] = -1;
1078             etCone+=etCell[ncell];
1079             nCellIn++;
1080           }
1081           else nCellOut++;
1082         }
1083
1084       // Select jets with et > background
1085       // estimate max fluctuation of background in cone
1086       Double_t ncellin = (Double_t)nCellIn;
1087       Double_t ntcell  = (Double_t)nCell;
1088       Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
1089       // min cone et
1090       Double_t etcmin = etCone ;  // could be used etCone - etmin !!
1091       //decisions !! etbmax < etcmin
1092
1093       for(Int_t mcell =0; mcell < nCell; mcell++){
1094         if(flagCell[mcell] == -1){
1095           if(etbmax < etcmin)
1096             flagCell[mcell] = 1; //flag cell as used
1097           else
1098             flagCell[mcell] = 0; // leave it free
1099         }
1100       }
1101       //store tmp jet info !!!
1102
1103       if(etbmax < etcmin) {
1104         etaAlgoJet[nJets]    = eta;
1105         phiAlgoJet[nJets]    = phi;
1106         etAlgoJet[nJets]     = etCone;
1107         ncellsAlgoJet[nJets] = nCellIn;
1108         nJets++;
1109       }
1110
1111     } // end of cells loop
1112
1113   //reorder jets by et in cone
1114   //sort jets by energy
1115   Int_t * idx  = new Int_t[nJets];
1116   TMath::Sort(nJets, etAlgoJet, idx);
1117   for(Int_t p = 0; p < nJets; p++)
1118     {
1119       etaJet[p] = etaAlgoJet[idx[p]];
1120       phiJet[p] = phiAlgoJet[idx[p]];
1121       etJet[p] = etAlgoJet[idx[p]];
1122       etallJet[p] = etAlgoJet[idx[p]];
1123       ncellsJet[p] = ncellsAlgoJet[idx[p]];
1124     }
1125
1126   //delete
1127   delete index;
1128   delete idx;
1129
1130 }
1131
1132 ////////////////////////////////////////////////////////////////////////
1133 void AliUA1JetFinderV2::SubtractBackg(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, const Float_t* ptT,
1134                                       const Int_t*vectT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
1135                                       const Float_t* cFlag2T, const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet,
1136                                       const Float_t* phiJet, Float_t* etsigJet, Int_t* multJet, Int_t* injet)
1137 {
1138   //
1139   // Background subtraction using cone method but without correction in dE/deta distribution
1140   // Cases to take into account the EMCal geometry are included
1141   //
1142
1143   //calculate energy inside and outside cones
1144   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1145   Int_t              fOpt = fReader->GetReaderHeader()->GetDetector();
1146   fDebug = fReader->GetReaderHeader()->GetDebug();
1147   Float_t rc= header->GetRadius();
1148   Float_t etIn[30];
1149   Float_t etOut = 0;
1150
1151   for(Int_t j=0;j<30;j++){etIn[j]=0.;}
1152
1153   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1154
1155     for(Int_t ijet=0; ijet<nJ; ijet++){
1156
1157       Float_t deta = etaT[jpart] - etaJet[ijet];
1158       Float_t dphi = phiT[jpart] - phiJet[ijet];
1159       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1160       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1161
1162       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1163       if(dr <= rc){ // particles inside this cone
1164         multJet[ijet]+=vectT[jpart];
1165         injet[jpart] = ijet;
1166
1167         if(cFlagT[jpart] == 1 || cFlag2T[jpart] == 1){ // pt cut
1168           etIn[ijet] += ptT[jpart];
1169           if(sFlagT[jpart] == 1) etsigJet[ijet]+= ptT[jpart];
1170         }
1171         break;
1172       }
1173     }// end jets loop
1174
1175     if(injet[jpart] == -1 && (cFlagT[jpart] == 1 || cFlag2T[jpart] == 1)){
1176       etOut += ptT[jpart]; // particle outside cones and pt cut
1177     }
1178   } //end particle loop
1179
1180   //estimate jets and background areas
1181   // TPC case
1182   if(fOpt == 0 || fOpt == 1){
1183     Float_t areaJet[30];
1184     Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1185
1186     for(Int_t k=0; k<nJ; k++){
1187       Float_t detamax = etaJet[k] + rc;
1188       Float_t detamin = etaJet[k] - rc;
1189       Float_t accmax = 0.0; Float_t accmin = 0.0;
1190       if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1191         Float_t h = header->GetLegoEtaMax() - etaJet[k];
1192         accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1193       }
1194       if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1195         Float_t h = header->GetLegoEtaMax() + etaJet[k];
1196         accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1197       }
1198       areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1199       areaOut = areaOut - areaJet[k];
1200     }
1201     //subtract background using area method
1202     for(Int_t ljet=0; ljet<nJ; ljet++){
1203       Float_t areaRatio = areaJet[ljet]/areaOut;
1204       etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1205     }
1206
1207     // estimate new total background
1208     Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1209     etbgTotalN = etOut*areaT/areaOut;
1210   }
1211   else { // If EMCal included
1212     Float_t areaJet[30];
1213     Float_t areaOut = 2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax() - header->GetLegoPhiMin());
1214     for(Int_t k=0; k<nJ; k++){
1215       Float_t detamax = etaJet[k] + rc;
1216       Float_t detamin = etaJet[k] - rc;
1217       Float_t dphimax = phiJet[k] + rc;
1218       Float_t dphimin = phiJet[k] - rc;
1219       Float_t eMax = header->GetLegoEtaMax();
1220       Float_t eMin = header->GetLegoEtaMin();
1221       Float_t pMax = header->GetLegoPhiMax();
1222       Float_t pMin = header->GetLegoPhiMin();
1223       Float_t accetamax = 0.0; Float_t accetamin = 0.0;
1224       Float_t accphimax = 0.0; Float_t accphimin = 0.0;
1225       if((detamax > eMax && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
1226          (detamax > eMax && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
1227         Float_t h = eMax - etaJet[k];
1228         accetamax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1229       }
1230       if((detamin < eMin && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
1231          (detamin < eMin && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
1232         Float_t h = eMax + etaJet[k];
1233         accetamin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1234       }
1235       if((dphimax > pMax && detamax >= (eMin+2*rc) && detamax <= eMax )||
1236          (dphimax > pMax && detamin <= (eMax-2*rc) && detamin >= eMin )){
1237         Float_t h = pMax - phiJet[k];
1238         accphimax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1239       }
1240       if((dphimin < eMin && detamax >= (eMin+2*rc) && detamax <= eMax )||
1241          (dphimin < eMin && detamin <= (eMax-2*rc) && detamin >= eMin )){
1242         Float_t h = phiJet[k] - pMin;
1243         accphimin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1244       }
1245
1246       if(detamax > eMax && dphimax > pMax ){
1247         Float_t he = eMax - etaJet[k];
1248         Float_t hp = pMax - phiJet[k];
1249         Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1250         Float_t alphae = TMath::ACos(he/rc);
1251         Float_t alphap = TMath::ACos(hp/rc);
1252         Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1253         if(rlim <= rc){
1254           accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1255           accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1256         }
1257         if(rlim > rc){
1258           accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1259           accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1260             ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1261             rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1262         }
1263       }
1264
1265       if(detamax > eMax && dphimin < pMin ){
1266         Float_t he = eMax - etaJet[k];
1267         Float_t hp = phiJet[k] - pMin;
1268         Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1269         Float_t alphae = TMath::ACos(he/rc);
1270         Float_t alphap = TMath::ACos(hp/rc);
1271         Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1272         if(rlim <= rc){
1273           accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1274           accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1275         }
1276         if(rlim > rc){
1277           accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1278           accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1279             ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1280             rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1281         }
1282       }
1283
1284       if(detamin < eMin && dphimax > pMax ){
1285         Float_t he = eMax + etaJet[k];
1286         Float_t hp = pMax - phiJet[k];
1287         Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1288         Float_t alphae = TMath::ACos(he/rc);
1289         Float_t alphap = TMath::ACos(hp/rc);
1290         Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1291         if(rlim <= rc){
1292           accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1293           accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1294         }
1295         if(rlim > rc){
1296           accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1297           accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1298             ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1299             rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1300         }
1301       }
1302
1303       if(detamin < eMin && dphimin < pMin ){
1304         Float_t he = eMax + etaJet[k];
1305         Float_t hp = phiJet[k] - pMin;
1306         Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1307         Float_t alphae = TMath::ACos(he/rc);
1308         Float_t alphap = TMath::ACos(hp/rc);
1309         Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1310         if(rlim <= rc){
1311           accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1312           accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1313         }
1314         if(rlim > rc){
1315           accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1316           accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1317             ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1318             rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1319         }
1320       }
1321       areaJet[k] = rc*rc*TMath::Pi() - accetamax - accetamin - accphimax - accphimin;
1322       areaOut = areaOut - areaJet[k];
1323     } // end loop on jets
1324
1325     //subtract background using area method
1326     for(Int_t ljet=0; ljet<nJ; ljet++){
1327       Float_t areaRatio = areaJet[ljet]/areaOut;
1328       etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1329     }
1330
1331     // estimate new total background
1332     Float_t areaT = 2*(header->GetLegoEtaMax()*header->GetLegoPhiMax());
1333     etbgTotalN = etOut*areaT/areaOut;
1334   }
1335
1336 }
1337
1338 ////////////////////////////////////////////////////////////////////////
1339 void AliUA1JetFinderV2::SubtractBackgC(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN,
1340                                        const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
1341                                        Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1342                                        Int_t* multJet, Int_t* injet)
1343 {
1344   //background subtraction using cone method but without correction in dE/deta distribution
1345
1346   //calculate energy inside and outside cones
1347   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1348   Float_t rc= header->GetRadius();
1349   Float_t etIn[30];
1350   Float_t etOut = 0;
1351   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1352     // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
1353     for(Int_t ijet=0; ijet<nJ; ijet++){
1354       Float_t deta = etaT[jpart] - etaJet[ijet];
1355       Float_t dphi = phiT[jpart] - phiJet[ijet];
1356       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1357       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1358       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1359       if(dr <= rc){ // particles inside this cone
1360         multJet[ijet]++;
1361         injet[jpart] = ijet;
1362         if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
1363           etIn[ijet] += ptT[jpart];
1364           if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
1365         }
1366         break;
1367       }
1368     }// end jets loop
1369     if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
1370       etOut += ptT[jpart]; // particle outside cones and pt cut
1371   } //end particle loop
1372
1373   //estimate jets and background areas
1374   Float_t areaJet[30];
1375   Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1376   for(Int_t k=0; k<nJ; k++){
1377     Float_t detamax = etaJet[k] + rc;
1378     Float_t detamin = etaJet[k] - rc;
1379     Float_t accmax = 0.0; Float_t accmin = 0.0;
1380     if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1381       Float_t h = header->GetLegoEtaMax() - etaJet[k];
1382       accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1383     }
1384     if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1385       Float_t h = header->GetLegoEtaMax() + etaJet[k];
1386       accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1387     }
1388     areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1389     areaOut = areaOut - areaJet[k];
1390   }
1391   //subtract background using area method
1392   for(Int_t ljet=0; ljet<nJ; ljet++){
1393     Float_t areaRatio = areaJet[ljet]/areaOut;
1394     etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1395   }
1396
1397   // estimate new total background
1398   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1399   etbgTotalN = etOut*areaT/areaOut;
1400
1401 }
1402
1403
1404 ////////////////////////////////////////////////////////////////////////
1405 void AliUA1JetFinderV2::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN,
1406                       const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
1407                       const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
1408                       Float_t* etsigJet, Int_t* multJet, Int_t* injet)
1409 {
1410
1411   //background subtraction using statistical method
1412   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1413   Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
1414
1415   //calculate energy inside
1416   Float_t rc= header->GetRadius();
1417   Float_t etIn[30];
1418
1419   for(Int_t jpart = 0; jpart < nIn; jpart++)
1420     { // loop for all particles in array
1421
1422       for(Int_t ijet=0; ijet<nJ; ijet++)
1423         {
1424           Float_t deta = etaT[jpart] - etaJet[ijet];
1425           Float_t dphi = phiT[jpart] - phiJet[ijet];
1426           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1427           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1428           Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1429           if(dr <= rc){ // particles inside this cone
1430             multJet[ijet]++;
1431             injet[jpart] = ijet;
1432             if(cFlagT[jpart] == 1){ // pt cut
1433               etIn[ijet]+= ptT[jpart];
1434               if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1435             }
1436             break;
1437           }
1438         }// end jets loop
1439     } //end particle loop
1440
1441   //calc jets areas
1442   Float_t areaJet[30];
1443   Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1444   for(Int_t k=0; k<nJ; k++)
1445     {
1446       Float_t detamax = etaJet[k] + rc;
1447       Float_t detamin = etaJet[k] - rc;
1448       Float_t accmax = 0.0; Float_t accmin = 0.0;
1449       if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1450         Float_t h = header->GetLegoEtaMax() - etaJet[k];
1451         accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1452       }
1453       if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1454       Float_t h = header->GetLegoEtaMax() + etaJet[k];
1455       accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1456       }
1457       areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1458     }
1459
1460   //subtract background using area method
1461   for(Int_t ljet=0; ljet<nJ; ljet++){
1462     Float_t areaRatio = areaJet[ljet]/areaOut;
1463     etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
1464   }
1465
1466   etbgTotalN = etbgStat;
1467 }
1468
1469 ////////////////////////////////////////////////////////////////////////
1470 void AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT,
1471                                           Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
1472                                           Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1473                                           Int_t* multJet, Int_t* injet)
1474 {
1475   // Cone background subtraction method taking into acount dEt/deta distribution
1476   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1477   //general
1478   Float_t rc= header->GetRadius();
1479   Float_t etamax = header->GetLegoEtaMax();
1480   Float_t etamin = header->GetLegoEtaMin();
1481   Int_t ndiv = 100;
1482
1483   // jet energy and area arrays
1484   TH1F* hEtJet[30];
1485   TH1F* hAreaJet[30];
1486   for(Int_t mjet=0; mjet<nJ; mjet++){
1487     char hEtname[256]; char hAreaname[256];
1488     sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
1489     hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
1490     hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
1491   }
1492   // background energy and area
1493   TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
1494   TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
1495
1496   //fill energies
1497   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1498     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
1499       Float_t deta = etaT[jpart] - etaJet[ijet];
1500       Float_t dphi = phiT[jpart] - phiJet[ijet];
1501       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1502       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1503       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1504       if(dr <= rc){ // particles inside this cone
1505         injet[jpart] = ijet;
1506         multJet[ijet]++;
1507         if(cFlagT[jpart] == 1){// pt cut
1508           hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
1509           if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1510         }
1511         break;
1512       }
1513     }// end jets loop
1514
1515     if(injet[jpart] == -1  && cFlagT[jpart] == 1)
1516       hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
1517   } //end particle loop
1518
1519   //calc areas
1520   Float_t eta0 = etamin;
1521   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
1522   Float_t eta1 = eta0 + etaw;
1523   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
1524     Float_t etac = eta0 + etaw/2.0;
1525     Float_t areabg = etaw*2.0*TMath::Pi();
1526     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
1527       Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
1528       Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
1529       Float_t acc0 = 0.0; Float_t acc1 = 0.0;
1530       Float_t areaj = 0.0;
1531       if(deta0 > rc && deta1 < rc){
1532         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1533         areaj = acc1;
1534       }
1535       if(deta0 < rc && deta1 > rc){
1536         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1537         areaj = acc0;
1538       }
1539       if(deta0 < rc && deta1 < rc){
1540         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1541         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1542         if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
1543         if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
1544         if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
1545       }
1546       hAreaJet[ijet]->Fill(etac,areaj);
1547       areabg = areabg - areaj;
1548     } // end jets loop
1549     hAreaBackg->Fill(etac,areabg);
1550     eta0 = eta1;
1551     eta1 = eta1 + etaw;
1552   } // end loop for all eta bins
1553
1554   //subtract background
1555   for(Int_t kjet=0; kjet<nJ; kjet++){
1556     etJet[kjet] = 0.0; // first  clear etJet for this jet
1557     for(Int_t bin = 0; bin< ndiv; bin++){
1558       if(hAreaJet[kjet]->GetBinContent(bin)){
1559         Float_t areab = hAreaBackg->GetBinContent(bin);
1560         Float_t etb = hEtBackg->GetBinContent(bin);
1561         Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
1562         etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
1563       }
1564     }
1565   }
1566
1567   // calc background total
1568   Double_t etOut = hEtBackg->Integral();
1569   Double_t areaOut = hAreaBackg->Integral();
1570   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1571   etbgTotalN = etOut*areaT/areaOut;
1572
1573   //delete
1574   for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
1575     delete hEtJet[ljet];
1576     delete hAreaJet[ljet];
1577   }
1578
1579   delete hEtBackg;
1580   delete hAreaBackg;
1581 }
1582
1583 ////////////////////////////////////////////////////////////////////////
1584 void AliUA1JetFinderV2::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
1585                       Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
1586                       Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1587                       Int_t* multJet, Int_t* injet)
1588 {
1589   // Ratio background subtraction method taking into acount dEt/deta distribution
1590   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1591   //factor F calc before
1592   Float_t bgRatioCut = header->GetBackgCutRatio();
1593
1594   //general
1595   Float_t rc= header->GetRadius();
1596   Float_t etamax = header->GetLegoEtaMax();
1597   Float_t etamin = header->GetLegoEtaMin();
1598   Int_t ndiv = 100;
1599
1600   // jet energy and area arrays
1601   TH1F* hEtJet[30];
1602   TH1F* hAreaJet[30];
1603   for(Int_t mjet=0; mjet<nJ; mjet++){
1604     char hEtname[256]; char hAreaname[256];
1605     sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
1606     hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);        // change range
1607     hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);  // change range
1608   }
1609   // background energy and area
1610   TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);         // change range
1611   TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);  // change range
1612
1613   //fill energies
1614   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1615     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
1616       Float_t deta = etaT[jpart] - etaJet[ijet];
1617       Float_t dphi = phiT[jpart] - phiJet[ijet];
1618       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1619       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1620       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1621       if(dr <= rc){ // particles inside this cone
1622         multJet[ijet]++;
1623         injet[jpart] = ijet;
1624         if(cFlagT[jpart] == 1){ //pt cut
1625           hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
1626           if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1627         }
1628         break;
1629       }
1630     }// end jets loop
1631     if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
1632   } //end particle loop
1633
1634   //calc areas
1635   Float_t eta0 = etamin;
1636   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
1637   Float_t eta1 = eta0 + etaw;
1638   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
1639     Float_t etac = eta0 + etaw/2.0;
1640     Float_t areabg = etaw*2.0*TMath::Pi();
1641     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
1642       Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
1643       Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
1644       Float_t acc0 = 0.0; Float_t acc1 = 0.0;
1645       Float_t areaj = 0.0;
1646       if(deta0 > rc && deta1 < rc){
1647         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1648         areaj = acc1;
1649       }
1650       if(deta0 < rc && deta1 > rc){
1651         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1652         areaj = acc0;
1653       }
1654       if(deta0 < rc && deta1 < rc){
1655         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1656         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1657         if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
1658         if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
1659         if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
1660       }
1661       hAreaJet[ijet]->Fill(etac,areaj);
1662       areabg = areabg - areaj;
1663     } // end jets loop
1664     hAreaBackg->Fill(etac,areabg);
1665     eta0 = eta1;
1666     eta1 = eta1 + etaw;
1667   } // end loop for all eta bins
1668
1669   //subtract background
1670   for(Int_t kjet=0; kjet<nJ; kjet++){
1671     etJet[kjet] = 0.0; // first  clear etJet for this jet
1672     for(Int_t bin = 0; bin< ndiv; bin++){
1673       if(hAreaJet[kjet]->GetBinContent(bin)){
1674         Float_t areab = hAreaBackg->GetBinContent(bin);
1675         Float_t etb = hEtBackg->GetBinContent(bin);
1676         Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
1677         etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
1678       }
1679     }
1680   }
1681
1682   // calc background total
1683   Double_t etOut = hEtBackg->Integral();
1684   Double_t areaOut = hAreaBackg->Integral();
1685   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1686   etbgTotalN = etOut*areaT/areaOut;
1687
1688   //delete
1689   for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
1690     delete hEtJet[ljet];
1691     delete hAreaJet[ljet];
1692   }
1693
1694   delete hEtBackg;
1695   delete hAreaBackg;
1696 }
1697
1698 ////////////////////////////////////////////////////////////////////////
1699 void AliUA1JetFinderV2::Reset()
1700 {
1701   fLego->Reset();
1702   fJets->ClearJets();
1703   AliJetFinder::Reset();
1704 }
1705
1706 ////////////////////////////////////////////////////////////////////////
1707 void AliUA1JetFinderV2::WriteJHeaderToFile()
1708 {
1709   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1710   header->Write();
1711 }
1712
1713 ////////////////////////////////////////////////////////////////////////
1714 void AliUA1JetFinderV2::InitTask(TChain* tree)
1715 {
1716
1717   // initializes some variables
1718   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1719   // book lego
1720   fLego = new TH2F("legoH","eta-phi",
1721                    header->GetLegoNbinEta(), header->GetLegoEtaMin(),
1722                    header->GetLegoEtaMax(),  header->GetLegoNbinPhi(),
1723                    header->GetLegoPhiMin(),  header->GetLegoPhiMax());
1724
1725   fDebug = fReader->GetReaderHeader()->GetDebug();
1726   fOpt = fReader->GetReaderHeader()->GetDetector();
1727
1728   // Tasks initialization
1729   if(fOpt>0)
1730     fReader->CreateTasks(tree);
1731
1732 }