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