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 #include "AliJetReaderHeader.h"
  38 #include "AliJetReader.h"
  39 #include "AliJetHeader.h"
  40
  41 class TArrayF;
  42 class TFile;
  43 class AliJetReader;
  44 class AliAODJet;
  45
  46 ClassImp(AliUA1JetFinderV2)
  47
  48
  49 ////////////////////////////////////////////////////////////////////////
  50 AliUA1JetFinderV2::AliUA1JetFinderV2() :
  51   AliJetFinder(),
  52   fLego(0),
  53   fDebug(0),
  54   fOpt(0)
  55 {
  56   //
  57   // Constructor
  58   //
  59 }
  60
  61 ////////////////////////////////////////////////////////////////////////
  62 AliUA1JetFinderV2::~AliUA1JetFinderV2()
  63 {
  64   //
  65   // Destructor
  66   //
  67 }
  68
  69 ////////////////////////////////////////////////////////////////////////
  70 void AliUA1JetFinderV2::FindJetsC()
  71 {
  72   //
  73   //  Used to find jets using charged particle momentum information
  74   //
  75   //  1) Fill cell map array
  76   //  2) calculate total energy and fluctuation level
  77   //  3) Run algorithm
  78   //     3.1) look centroides in cell map
  79   //     3.2) calculate total energy in cones
  80   //     3.3) flag as a possible jet
  81   //     3.4) reorder cones by energy
  82   //  4) subtract backg in accepted jets
  83   //  5) fill AliJet list
  84
  85   //  Transform input to pt,eta,phi plus lego
  86
  87   AliUA1JetHeaderV1* header  = (AliUA1JetHeaderV1*) fHeader;
  88   TClonesArray*      lvArray = fReader->GetMomentumArray();
  89   Int_t              nIn     = lvArray->GetEntries();
  90
  91   if (nIn == 0) return;
  92
  93   // local arrays for input
  94   Float_t* ptT    = new Float_t[nIn];
  95   Float_t* etaT   = new Float_t[nIn];
  96   Float_t* phiT   = new Float_t[nIn];
  97   Float_t* cFlagT = new Float_t[nIn]; // Temporarily added
  98   Float_t* sFlagT = new Float_t[nIn]; // Temporarily added
  99   Int_t*   injet  = new Int_t[nIn];
 100
 101   //total energy in array
 102   Float_t  etbgTotal = 0.0;
 103   TH1F*    hPtTotal  = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
 104
 105   // load input vectors and calculate total energy in array
 106   for (Int_t i = 0; i < nIn; i++){
 107     TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
 108     ptT[i]  = lv->Pt();
 109     etaT[i] = lv->Eta();
 110     phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
 111     cFlagT[i] = fReader->GetCutFlag(i);
 112     sFlagT[i] = fReader->GetSignalFlag(i);
 113
 114     if (fReader->GetCutFlag(i) != 1) continue;
 115     fLego->Fill(etaT[i], phiT[i], ptT[i]);
 116     hPtTotal->Fill(ptT[i]);
 117     etbgTotal+= ptT[i];
 118   }
 119
 120
 121   // calculate total energy and fluctuation in map
 122   Double_t meanpt   = hPtTotal->GetMean();
 123   Double_t ptRMS    = hPtTotal->GetRMS();
 124   Double_t npart    = hPtTotal->GetEntries();
 125   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
 126
 127   // arrays to hold jets
 128   Float_t* etaJet    = new Float_t[30];  // eta jet
 129   Float_t* phiJet    = new Float_t[30];  // phi jet
 130   Float_t* etJet     = new Float_t[30];  // et jet
 131   Float_t* etsigJet  = new Float_t[30];  // signal et in jet
 132   Float_t* etallJet  = new Float_t[30];  // total et in jet (tmp variable)
 133   Int_t*   ncellsJet = new Int_t[30];
 134   Int_t*   multJet   = new Int_t[30];
 135   //--- Added for jet reordering at the end of the jet finding procedure
 136   Float_t* etaJetOk    = new Float_t[30];
 137   Float_t* phiJetOk    = new Float_t[30];
 138   Float_t* etJetOk     = new Float_t[30];
 139   Float_t* etsigJetOk  = new Float_t[30];  // signal et in jet
 140   Float_t* etallJetOk  = new Float_t[30];  // total et in jet (tmp variable)
 141   Int_t*   ncellsJetOk = new Int_t[30];
 142   Int_t*   multJetOk   = new Int_t[30];
 143   //--------------------------
 144   Int_t nJets; // to hold number of jets found by algorithm
 145   Int_t nj;    // number of jets accepted
 146   Float_t prec  = header->GetPrecBg();
 147   Float_t bgprec = 1;
 148   while(bgprec > prec){
 149     //reset jet arrays in memory
 150     memset(etaJet,0,sizeof(Float_t)*30);
 151     memset(phiJet,0,sizeof(Float_t)*30);
 152     memset(etJet,0,sizeof(Float_t)*30);
 153     memset(etallJet,0,sizeof(Float_t)*30);
 154     memset(etsigJet,0,sizeof(Float_t)*30);
 155     memset(ncellsJet,0,sizeof(Int_t)*30);
 156     memset(multJet,0,sizeof(Int_t)*30);
 157     //--- Added for jet reordering at the end of the jet finding procedure
 158     memset(etaJetOk,0,sizeof(Float_t)*30);
 159     memset(phiJetOk,0,sizeof(Float_t)*30);
 160     memset(etJetOk,0,sizeof(Float_t)*30);
 161     memset(etallJetOk,0,sizeof(Float_t)*30);
 162     memset(etsigJetOk,0,sizeof(Float_t)*30);
 163     memset(ncellsJetOk,0,sizeof(Int_t)*30);
 164     memset(multJetOk,0,sizeof(Int_t)*30);
 165     //--------------------------
 166     nJets = 0;
 167     nj = 0;
 168
 169     // reset particles-jet array in memory
 170     memset(injet,-1,sizeof(Int_t)*nIn);
 171     //run cone algorithm finder
 172     RunAlgoritmC(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
 173
 174     //run background subtraction
 175     if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
 176       nj = header->GetNAcceptJets();
 177     else
 178       nj = nJets;
 179     //subtract background
 180     Float_t etbgTotalN = 0.0; //new background
 181     if(header->GetBackgMode() == 1) // standard
 182       SubtractBackgC(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 183     if(header->GetBackgMode() == 2) //cone
 184       SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 185     if(header->GetBackgMode() == 3) //ratio
 186       SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 187     if(header->GetBackgMode() == 4) //statistic
 188       SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 189     //calc precision
 190     if(etbgTotalN != 0.0)
 191       bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
 192     else
 193       bgprec = 0;
 194     etbgTotal = etbgTotalN; // update with new background estimation
 195   } //end while
 196
 197   // add jets to list
 198   Int_t* idxjets = new Int_t[nj];
 199   Int_t nselectj = 0;
 200   printf("Found %d jets \n", nj);
 201
 202   // Reorder jets by et in cone
 203   Int_t * idx  = new Int_t[nJets];
 204   TMath::Sort(nJets, etJet, idx);
 205   for(Int_t p = 0; p < nJets; p++){
 206     etaJetOk[p]    = etaJet[idx[p]];
 207     phiJetOk[p]    = phiJet[idx[p]];
 208     etJetOk[p]     = etJet[idx[p]];
 209     etallJetOk[p]  = etJet[idx[p]];
 210     etsigJetOk[p]  = etsigJet[idx[p]];
 211     ncellsJetOk[p] = ncellsJet[idx[p]];
 212     multJetOk[p]   = multJet[idx[p]];
 213   }
 214
 215   for(Int_t kj=0; kj<nj; kj++)
 216     {
 217       if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
 218           (etaJetOk[kj] < (header->GetJetEtaMin())) ||
 219           (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
 220       Float_t px, py,pz,en; // convert to 4-vector
 221       px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
 222       py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
 223       pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
 224       en = TMath::Sqrt(px * px + py * py + pz * pz);
 225
 226       AliAODJet jet(px, py, pz, en);
 227       jet.Print("");
 228
 229       AddJet(jet);
 230
 231       idxjets[nselectj] = kj;
 232       nselectj++;
 233     }
 234
 235   //add signal percentage and total signal  in AliJets for analysis tool
 236   Float_t* percentage  = new Float_t[nselectj];
 237   Int_t* ncells      = new Int_t[nselectj];
 238   Int_t* mult        = new Int_t[nselectj];
 239   for(Int_t i = 0; i< nselectj; i++)
 240     {
 241       percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
 242       ncells[i] = ncellsJetOk[idxjets[i]];
 243       mult[i] = multJetOk[idxjets[i]];
 244     }
 245
 246   //add particle-injet relationship ///
 247   for(Int_t bj = 0; bj < nIn; bj++)
 248     {
 249       if(injet[bj] == -1) continue; //background particle
 250       Int_t bflag = 0;
 251       for(Int_t ci = 0; ci< nselectj; ci++){
 252         if(injet[bj] == idxjets[ci]){
 253           injet[bj]= ci;
 254           bflag++;
 255           break;
 256         }
 257       }
 258       if(bflag == 0) injet[bj] = -1; // set as background particle
 259     }
 260
 261
 262   //delete
 263   delete[] ptT;
 264   delete[] etaT;
 265   delete[] phiT;
 266   delete[] cFlagT;
 267   delete[] sFlagT;
 268   delete[] injet;
 269   delete[] hPtTotal;
 270   delete[] etaJet;
 271   delete[] phiJet;
 272   delete[] etJet;
 273   delete[] etsigJet;
 274   delete[] etallJet;
 275   delete[] ncellsJet;
 276   delete[] multJet;
 277   delete[] idxjets;
 278   delete[] percentage;
 279   delete[] ncells;
 280   delete[] mult;
 281   //--- Added for jet reordering
 282   delete etaJetOk;
 283   delete phiJetOk;
 284   delete etJetOk;
 285   delete etsigJetOk;
 286   delete etallJetOk;
 287   delete ncellsJetOk;
 288   delete multJetOk;
 289   //--------------------------
 290
 291 }
 292
 293 ////////////////////////////////////////////////////////////////////////
 294 void AliUA1JetFinderV2::FindJets()
 295 {
 296   //
 297   //  Used to find jets using charged particle momentum information
 298   //  & neutral energy from calo cells
 299   //
 300   //  1) Fill cell map array
 301   //  2) calculate total energy and fluctuation level
 302   //  3) Run algorithm
 303   //     3.1) look centroides in cell map
 304   //     3.2) calculate total energy in cones
 305   //     3.3) flag as a possible jet
 306   //     3.4) reorder cones by energy
 307   //  4) subtract backg in accepted jets
 308   //  5) fill AliJet list
 309
 310   // transform input to pt,eta,phi plus lego
 311
 312   AliUA1JetHeaderV1* header   = (AliUA1JetHeaderV1*) fHeader;
 313   TClonesArray*      fUnit    = fReader->GetUnitArray();
 314   Int_t              nCand    = fReader->GetNumCandidate();
 315   Int_t              nCandCut = fReader->GetNumCandidateCut();
 316   Int_t              nIn      = fUnit->GetEntries();
 317   Float_t            ptMin   = fReader->GetReaderHeader()->GetPtCut();
 318
 319   fDebug   = fReader->GetReaderHeader()->GetDebug();
 320   if (nIn == 0) return;
 321
 322   Int_t nCandidateCut = 0;
 323   Int_t nCandidate = 0;
 324
 325   nCandidate = nCand;
 326   nCandidateCut = nCandCut;
 327
 328   // local arrays for input No Cuts
 329   // Both pt < ptMin and pt > ptMin
 330   Float_t*   ptT       = new Float_t[nCandidate];
 331   Float_t*   en2T      = new Float_t[nCandidate];
 332   Float_t*   pt2T      = new Float_t[nCandidate];
 333   Int_t*     detT      = new Int_t[nCandidate];
 334   Float_t*   etaT      = new Float_t[nCandidate];
 335   Float_t*   phiT      = new Float_t[nCandidate];
 336   Float_t*   cFlagT    = new Float_t[nCandidate];
 337   Float_t*   cFlag2T   = new Float_t[nCandidate];
 338   Float_t*   sFlagT    = new Float_t[nCandidate];
 339   Float_t*   cClusterT = new Float_t[nCandidate];
 340   Int_t*     vectT     = new Int_t[nCandidate];
 341   Int_t      loop1     = 0;
 342   Int_t*     injet     = new Int_t[nCandidate];
 343   Int_t*     sflag     = new Int_t[nCandidate];
 344   TRefArray* trackRef  = new TRefArray();
 345
 346   //total energy in array
 347   Float_t  etbgTotal = 0.0;
 348   TH1F*    hPtTotal  = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
 349
 350   // Input cell info
 351   Float_t *etCell    = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
 352   Float_t *etaCell   = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
 353   Float_t *phiCell   = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
 354   Int_t   *flagCell  = new Int_t[nIn];   //! Cell phi - Extracted from UnitArray
 355   Float_t *etCell2   = new Float_t[nIn];  //! Cell Energy - Extracted from UnitArray
 356   Float_t *etaCell2  = new Float_t[nIn]; //! Cell eta - Extracted from UnitArray
 357   Float_t *phiCell2  = new Float_t[nIn]; //! Cell phi - Extracted from UnitArray
 358   Int_t   *flagCell2 = new Int_t[nIn];  //! Cell phi - Extracted from UnitArray
 359
 360   // Information extracted from fUnitArray
 361   // Load input vectors and calculate total energy in array
 362   for(Int_t i=0; i<nIn; i++)
 363     {
 364       // Recover particle information from UnitArray
 365
 366       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(i);
 367       TRefArray* ref = uArray->GetUnitTrackRef();
 368       Int_t nRef = ref->GetEntries();
 369
 370       if(uArray->GetUnitEnergy()>0.){
 371
 372         for(Int_t jpart=0; jpart<nRef;jpart++)
 373           trackRef->Add((AliVTrack*)ref->At(jpart));
 374         ptT[loop1]   = uArray->GetUnitEnergy();
 375         detT[loop1]  = uArray->GetUnitDetectorFlag();
 376         etaT[loop1]  = uArray->GetUnitEta();
 377         phiT[loop1]  = uArray->GetUnitPhi();
 378         cFlagT[loop1]= uArray->GetUnitCutFlag();   // pt cut tpc
 379         cFlag2T[loop1]= uArray->GetUnitCutFlag2(); // pt cut emcal
 380         sFlagT[loop1]= uArray->GetUnitSignalFlag();
 381         vectT[loop1] = nRef;
 382         if(cFlagT[loop1] == 1 || cFlag2T[loop1] == 1) {
 383           pt2T[loop1] = 0.;
 384           en2T[loop1] = 0.;
 385           if(detT[loop1]==1){
 386             en2T[loop1] = ptT[loop1] - header->GetMinCellEt();
 387             if(en2T[loop1] < 0) en2T[loop1]=0;
 388             hPtTotal->Fill(en2T[loop1]);
 389             etbgTotal += en2T[loop1];
 390           }
 391           if(detT[loop1]==0){ // TPC+ITS
 392             Float_t pt = 0.;
 393             for(Int_t j=0; j<nRef;j++){
 394               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
 395               x = ((AliVTrack*)ref->At(j))->Px();
 396               y = ((AliVTrack*)ref->At(j))->Py();
 397               z = ((AliVTrack*)ref->At(j))->Pz();
 398               pt = TMath::Sqrt(x*x+y*y);
 399               if(pt>ptMin) {
 400                 pt2T[loop1] += pt;
 401                 en2T[loop1] += pt;
 402                 hPtTotal->Fill(pt);
 403                 etbgTotal+= pt;
 404               }
 405             }
 406           }
 407           if(detT[loop1]==2) { // EMCal
 408             Float_t ptCTot = 0.;
 409             Float_t pt = 0.;
 410             Float_t enC = 0.;
 411             for(Int_t j=0; j<nRef;j++){
 412               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
 413               x = ((AliVTrack*)ref->At(j))->Px();
 414               y = ((AliVTrack*)ref->At(j))->Py();
 415               z = ((AliVTrack*)ref->At(j))->Pz();
 416               pt = TMath::Sqrt(x*x+y*y);
 417               if(pt>ptMin) {
 418                 pt2T[loop1]+=pt;
 419                 en2T[loop1]+=pt;
 420                 hPtTotal->Fill(pt);
 421                 etbgTotal+= pt;
 422               }
 423               ptCTot += pt;
 424             }
 425             enC = ptT[loop1] - ptCTot - header->GetMinCellEt();
 426             if(enC < 0.) enC=0.;
 427             en2T[loop1] += enC;
 428             hPtTotal->Fill(enC);
 429             etbgTotal+= enC;
 430           }
 431         }
 432         loop1++;
 433       }
 434
 435       if(uArray->GetUnitCutFlag()==1) {
 436         if(uArray->GetUnitDetectorFlag()==1){ // EMCal case
 437           etCell[i] = uArray->GetUnitEnergy() - header->GetMinCellEt();
 438           if ((uArray->GetUnitEnergy() - header->GetMinCellEt()) < 0.0) etCell[i]=0.;
 439           etaCell[i] = uArray->GetUnitEta();
 440           phiCell[i] = uArray->GetUnitPhi();
 441           flagCell[i] = 0; // default
 442           etCell2[i] = etCell[i];
 443           etaCell2[i] = uArray->GetUnitEta();
 444           phiCell2[i] = uArray->GetUnitPhi();
 445           flagCell2[i] = 0; // default
 446         }
 447         if(uArray->GetUnitDetectorFlag()==0){ // TPC case
 448           Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
 449           for(Int_t j=0; j<nRef;j++)
 450             {
 451               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
 452               x = ((AliVTrack*)ref->At(j))->Px();
 453               y = ((AliVTrack*)ref->At(j))->Py();
 454               z = ((AliVTrack*)ref->At(j))->Pz();
 455               pt = TMath::Sqrt(x*x+y*y);
 456               if(pt>ptMin) {
 457                 et1 += pt;
 458                 et2 += pt;
 459               }
 460             }
 461           etCell[i] = et1;
 462           etCell2[i] = et2;
 463           if(et1 < 0.) etCell[i] = etCell2[i] = 0.;
 464           etaCell[i] = uArray->GetUnitEta();
 465           phiCell[i] = uArray->GetUnitPhi();
 466           flagCell[i] = 0; // default
 467           etaCell2[i] = uArray->GetUnitEta();
 468           phiCell2[i] = uArray->GetUnitPhi();
 469           flagCell2[i] = 0; // default
 470         }
 471         if(uArray->GetUnitDetectorFlag()==2){ // TPC + EMCal case
 472           Float_t ptCTot = 0.;
 473           Float_t pt = 0.; Float_t et1 = 0.; Float_t et2 = 0.;
 474           Float_t enC = 0.;
 475           for(Int_t j=0; j<nRef;j++)
 476             {
 477               Float_t x=0.;  Float_t y=0.;  Float_t z=0.;
 478               x = ((AliVTrack*)ref->At(j))->Px();
 479               y = ((AliVTrack*)ref->At(j))->Py();
 480               z = ((AliVTrack*)ref->At(j))->Pz();
 481               pt = TMath::Sqrt(x*x+y*y);
 482               if(pt>ptMin) {
 483                 et1 += pt;
 484                 et2 += pt;
 485               }
 486               ptCTot += pt;
 487             }
 488           enC = uArray->GetUnitEnergy() - ptCTot;
 489           etCell[i] = et1 + enC - header->GetMinCellEt();
 490           etCell2[i] = et2 + enC - header->GetMinCellEt();
 491           if((enC + et1 - header->GetMinCellEt()) < 0.) etCell[i] = etCell2[i] = 0.;
 492           etaCell[i] = uArray->GetUnitEta();
 493           phiCell[i] = uArray->GetUnitPhi();
 494           flagCell[i] = 0; // default
 495           etaCell2[i] = uArray->GetUnitEta();
 496           phiCell2[i] = uArray->GetUnitPhi();
 497           flagCell2[i] = 0; // default
 498         }
 499       }
 500       else {
 501         etCell[i]  = 0.;
 502         etaCell[i] = uArray->GetUnitEta();
 503         phiCell[i] = uArray->GetUnitPhi();
 504         flagCell[i] = 0;
 505         etCell2[i]  = 0.;
 506         etaCell2[i] = uArray->GetUnitEta();
 507         phiCell2[i] = uArray->GetUnitPhi();
 508         flagCell2[i] = 0;
 509       }
 510     } // end loop on nCandidate
 511
 512
 513   // calculate total energy and fluctuation in map
 514   Double_t meanpt = hPtTotal->GetMean();
 515   Double_t ptRMS = hPtTotal->GetRMS();
 516   Double_t npart = hPtTotal->GetEntries();
 517   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
 518
 519   // arrays to hold jets
 520   Float_t* etaJet    = new Float_t[30];
 521   Float_t* phiJet    = new Float_t[30];
 522   Float_t* etJet     = new Float_t[30];
 523   Float_t* etsigJet  = new Float_t[30]; //signal et in jet
 524   Float_t* etallJet  = new Float_t[30];  // total et in jet (tmp variable)
 525   Int_t*   ncellsJet = new Int_t[30];
 526   Int_t*   multJet   = new Int_t[30];
 527   //--- Added by me for jet reordering at the end of the jet finding procedure
 528   Float_t* etaJetOk    = new Float_t[30];
 529   Float_t* phiJetOk    = new Float_t[30];
 530   Float_t* etJetOk     = new Float_t[30];
 531   Float_t* etsigJetOk  = new Float_t[30]; //signal et in jet
 532   Float_t* etallJetOk  = new Float_t[30];  // total et in jet (tmp variable)
 533   Int_t*   ncellsJetOk = new Int_t[30];
 534   Int_t*   multJetOk   = new Int_t[30];
 535   //--------------------------
 536   Int_t    nJets; // to hold number of jets found by algorithm
 537   Int_t    nj;    // number of jets accepted
 538   Float_t  prec  = header->GetPrecBg();
 539   Float_t  bgprec = 1;
 540
 541   while(bgprec > prec){
 542
 543     //reset jet arrays in memory
 544     memset(etaJet,0,sizeof(Float_t)*30);
 545     memset(phiJet,0,sizeof(Float_t)*30);
 546     memset(etJet,0,sizeof(Float_t)*30);
 547     memset(etallJet,0,sizeof(Float_t)*30);
 548     memset(etsigJet,0,sizeof(Float_t)*30);
 549     memset(ncellsJet,0,sizeof(Int_t)*30);
 550     memset(multJet,0,sizeof(Int_t)*30);
 551     //--- Added by me for jet reordering at the end of the jet finding procedure
 552     memset(etaJetOk,0,sizeof(Float_t)*30);
 553     memset(phiJetOk,0,sizeof(Float_t)*30);
 554     memset(etJetOk,0,sizeof(Float_t)*30);
 555     memset(etallJetOk,0,sizeof(Float_t)*30);
 556     memset(etsigJetOk,0,sizeof(Float_t)*30);
 557     memset(ncellsJetOk,0,sizeof(Int_t)*30);
 558     memset(multJetOk,0,sizeof(Int_t)*30);
 559
 560     nJets = 0;
 561     nj = 0;
 562
 563     // reset particles-jet array in memory
 564     memset(injet,-1,sizeof(Int_t)*nCandidate);
 565     //run cone algorithm finder
 566     RunAlgoritm(nIn,etCell,etaCell,phiCell,flagCell,etCell2,etaCell2,phiCell2,
 567                 flagCell2,etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,
 568                 etallJet,ncellsJet);
 569
 570     //run background subtraction
 571     if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
 572       nj = header->GetNAcceptJets();
 573     else
 574       nj = nJets;
 575
 576     //subtract background
 577     Float_t etbgTotalN = 0.0; //new background
 578     if(header->GetBackgMode() == 1) // standard
 579       SubtractBackg(nCandidate,nj,etbgTotalN,en2T,vectT,etaT,phiT,cFlagT,cFlag2T,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 580     // To be modified ------------------------
 581     if(header->GetBackgMode() == 2) //cone
 582       SubtractBackgCone(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 583     if(header->GetBackgMode() == 3) //ratio
 584       SubtractBackgRatio(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 585     if(header->GetBackgMode() == 4) //statistic
 586       SubtractBackgStat(nCandidate,nj,etbgTotalN,ptT,etaT,phiT,cFlagT,sFlagT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
 587     //----------------------------------------
 588     //calc precision
 589     if(etbgTotalN != 0.0)
 590       bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
 591     else
 592       bgprec = 0;
 593     etbgTotal = etbgTotalN; // update with new background estimation
 594   } //end while
 595
 596   // add jets to list
 597   Int_t* idxjets = new Int_t[nj];
 598   Int_t nselectj = 0;
 599   printf("Found %d jets \n", nj);
 600
 601   // Reorder jets by et in cone
 602   // Sort jets by energy
 603   Int_t * idx  = new Int_t[nJets];
 604   TMath::Sort(nJets, etJet, idx);
 605   for(Int_t p = 0; p < nJets; p++)
 606     {
 607       etaJetOk[p]    = etaJet[idx[p]];
 608       phiJetOk[p]    = phiJet[idx[p]];
 609       etJetOk[p]     = etJet[idx[p]];
 610       etallJetOk[p]  = etJet[idx[p]];
 611       etsigJetOk[p]  = etsigJet[idx[p]];
 612       ncellsJetOk[p] = ncellsJet[idx[p]];
 613       multJetOk[p]   = multJet[idx[p]];
 614     }
 615
 616   TRefArray *refs = 0;
 617   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
 618   if (fromAod) refs = fReader->GetReferences();
 619   Int_t nTracks = 0;
 620   if (fromAod) nTracks = ((TRefArray*)refs)->GetEntries();
 621   Int_t* trackinjet     = new Int_t[nTracks];
 622   for(Int_t it=0; it<nTracks; it++) trackinjet[it]=-1;
 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
 635       AliAODJet jet(px, py, pz, en);
 636       jet.Print("");
 637
 638       if (fromAod){
 639         for(Int_t jpart = 0; jpart < nTracks; jpart++) { // loop for all particles in array
 640           Float_t deta = ((AliAODTrack*)refs->At(jpart))->Eta() - etaJetOk[kj];
 641           Float_t dphi = ((AliAODTrack*)refs->At(jpart))->Phi() - phiJetOk[kj];
 642           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 643           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 644
 645           Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 646           if(dr <= header->GetRadius() && fReader->GetCutFlag(jpart) == 1) // particles inside this cone
 647             if(trackinjet[jpart]==-1) trackinjet[jpart] = kj;
 648             else if(fDebug>10) printf("The track already belongs to jet %d \n",trackinjet[jpart]);
 649           if(trackinjet[jpart]==kj)
 650             jet.AddTrack(refs->At(jpart));  // check if the particle belongs to the jet and add the ref
 651         }
 652       }
 653
 654       AddJet(jet);
 655
 656       idxjets[nselectj] = kj;
 657       nselectj++;
 658     }
 659
 660   //add signal percentage and total signal  in AliJets for analysis tool
 661   Float_t* percentage  = new Float_t[nselectj];
 662   Int_t* ncells      = new Int_t[nselectj];
 663   Int_t* mult        = new Int_t[nselectj];
 664   for(Int_t i = 0; i< nselectj; i++)
 665     {
 666       percentage[i] = etsigJetOk[idxjets[i]]/etJetOk[idxjets[i]];
 667       ncells[i]     = ncellsJetOk[idxjets[i]];
 668       mult[i]       = multJetOk[idxjets[i]];
 669     }
 670
 671   //add particle-injet relationship ///
 672   for(Int_t bj = 0; bj < nCandidate; bj++)
 673     {
 674       if(injet[bj] == -1) continue; //background particle
 675       Int_t bflag = 0;
 676       for(Int_t ci = 0; ci< nselectj; ci++){
 677         if(injet[bj] == idxjets[ci]){
 678           injet[bj]= ci;
 679           bflag++;
 680           break;
 681         }
 682       }
 683     if(bflag == 0) injet[bj] = -1; // set as background particle
 684   }
 685
 686
 687   //delete
 688   delete ptT;
 689   delete en2T;
 690   delete pt2T;
 691   delete etaT;
 692   delete phiT;
 693   trackRef->Delete();
 694   delete trackRef;
 695   delete detT;
 696   delete cFlagT;
 697   delete cFlag2T;
 698   delete sFlagT;
 699   delete cClusterT;
 700   delete vectT;
 701   delete injet;
 702   delete sflag;
 703   delete hPtTotal;
 704   delete etCell;
 705   delete etaCell;
 706   delete phiCell;
 707   delete flagCell;
 708   delete etCell2;
 709   delete etaCell2;
 710   delete phiCell2;
 711   delete flagCell2;
 712   delete etaJet;
 713   delete phiJet;
 714   delete etJet;
 715   delete etsigJet;
 716   delete etallJet;
 717   delete ncellsJet;
 718   delete multJet;
 719   //--- Added for jet reordering
 720   delete etaJetOk;
 721   delete phiJetOk;
 722   delete etJetOk;
 723   delete etsigJetOk;
 724   delete etallJetOk;
 725   delete ncellsJetOk;
 726   delete multJetOk;
 727   //--------------------------
 728   delete trackinjet;
 729   delete idxjets;
 730   delete percentage;
 731   delete ncells;
 732   delete mult;
 733
 734 }
 735
 736 ////////////////////////////////////////////////////////////////////////
 737 void AliUA1JetFinderV2::RunAlgoritm(Int_t nIn, Float_t* etCell, Float_t* etaCell, Float_t* phiCell,
 738                                     Int_t* flagCell, const Float_t* etCell2, const Float_t* etaCell2,
 739                                     const Float_t* phiCell2, const Int_t* flagCell2, Float_t etbgTotal,
 740                                     Double_t dEtTotal, Int_t& nJets, Float_t* etJet,Float_t* etaJet,
 741                                     Float_t* phiJet, Float_t* etallJet, Int_t* ncellsJet)
 742 {
 743   //
 744   // Main method for jet finding
 745   // UA1 base cone finder
 746   //
 747
 748   Int_t nCell  = nIn;
 749   fDebug = fReader->GetReaderHeader()->GetDebug();
 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(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN, const Float_t* ptT,
1141                                       const Int_t*vectT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
1142                                       const Float_t* cFlag2T, const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet,
1143                                       const Float_t* phiJet, 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   Int_t              fOpt = fReader->GetReaderHeader()->GetDetector();
1153   fDebug = fReader->GetReaderHeader()->GetDebug();
1154   Float_t rc= header->GetRadius();
1155   Float_t etIn[30];
1156   Float_t etOut = 0;
1157
1158   for(Int_t j=0;j<30;j++){etIn[j]=0.;}
1159
1160   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1161
1162     for(Int_t ijet=0; ijet<nJ; ijet++){
1163
1164       Float_t deta = etaT[jpart] - etaJet[ijet];
1165       Float_t dphi = phiT[jpart] - phiJet[ijet];
1166       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1167       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1168
1169       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1170       if(dr <= rc){ // particles inside this cone
1171         multJet[ijet]+=vectT[jpart];
1172         injet[jpart] = ijet;
1173
1174         if(cFlagT[jpart] == 1 || cFlag2T[jpart] == 1){ // pt cut
1175           etIn[ijet] += ptT[jpart];
1176           if(sFlagT[jpart] == 1) etsigJet[ijet]+= ptT[jpart];
1177         }
1178         break;
1179       }
1180     }// end jets loop
1181
1182     if(injet[jpart] == -1 && (cFlagT[jpart] == 1 || cFlag2T[jpart] == 1)){
1183       etOut += ptT[jpart]; // particle outside cones and pt cut
1184     }
1185   } //end particle loop
1186
1187   //estimate jets and background areas
1188   // TPC case
1189   if(fOpt == 0 || fOpt == 1){
1190     Float_t areaJet[30];
1191     Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1192
1193     for(Int_t k=0; k<nJ; k++){
1194       Float_t detamax = etaJet[k] + rc;
1195       Float_t detamin = etaJet[k] - rc;
1196       Float_t accmax = 0.0; Float_t accmin = 0.0;
1197       if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1198         Float_t h = header->GetLegoEtaMax() - etaJet[k];
1199         accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1200       }
1201       if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1202         Float_t h = header->GetLegoEtaMax() + etaJet[k];
1203         accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1204       }
1205       areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1206       areaOut = areaOut - areaJet[k];
1207     }
1208     //subtract background using area method
1209     for(Int_t ljet=0; ljet<nJ; ljet++){
1210       Float_t areaRatio = areaJet[ljet]/areaOut;
1211       etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1212     }
1213
1214     // estimate new total background
1215     Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1216     etbgTotalN = etOut*areaT/areaOut;
1217   }
1218   else { // If EMCal included
1219     Float_t areaJet[30];
1220     Float_t areaOut = 2*(header->GetLegoEtaMax())*(header->GetLegoPhiMax() - header->GetLegoPhiMin());
1221     for(Int_t k=0; k<nJ; k++){
1222       Float_t detamax = etaJet[k] + rc;
1223       Float_t detamin = etaJet[k] - rc;
1224       Float_t dphimax = phiJet[k] + rc;
1225       Float_t dphimin = phiJet[k] - rc;
1226       Float_t eMax = header->GetLegoEtaMax();
1227       Float_t eMin = header->GetLegoEtaMin();
1228       Float_t pMax = header->GetLegoPhiMax();
1229       Float_t pMin = header->GetLegoPhiMin();
1230       Float_t accetamax = 0.0; Float_t accetamin = 0.0;
1231       Float_t accphimax = 0.0; Float_t accphimin = 0.0;
1232       if((detamax > eMax && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
1233          (detamax > eMax && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
1234         Float_t h = eMax - etaJet[k];
1235         accetamax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1236       }
1237       if((detamin < eMin && dphimax >= (pMin+2*rc) && dphimax <= pMax )||
1238          (detamin < eMin && dphimin <= (pMax-2*rc) && dphimin >= pMin )){
1239         Float_t h = eMax + etaJet[k];
1240         accetamin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1241       }
1242       if((dphimax > pMax && detamax >= (eMin+2*rc) && detamax <= eMax )||
1243          (dphimax > pMax && detamin <= (eMax-2*rc) && detamin >= eMin )){
1244         Float_t h = pMax - phiJet[k];
1245         accphimax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1246       }
1247       if((dphimin < eMin && detamax >= (eMin+2*rc) && detamax <= eMax )||
1248          (dphimin < eMin && detamin <= (eMax-2*rc) && detamin >= eMin )){
1249         Float_t h = phiJet[k] - pMin;
1250         accphimin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1251       }
1252
1253       if(detamax > eMax && dphimax > pMax ){
1254         Float_t he = eMax - etaJet[k];
1255         Float_t hp = pMax - phiJet[k];
1256         Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1257         Float_t alphae = TMath::ACos(he/rc);
1258         Float_t alphap = TMath::ACos(hp/rc);
1259         Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1260         if(rlim <= rc){
1261           accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1262           accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1263         }
1264         if(rlim > rc){
1265           accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1266           accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1267             ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1268             rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1269         }
1270       }
1271
1272       if(detamax > eMax && dphimin < pMin ){
1273         Float_t he = eMax - etaJet[k];
1274         Float_t hp = phiJet[k] - pMin;
1275         Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1276         Float_t alphae = TMath::ACos(he/rc);
1277         Float_t alphap = TMath::ACos(hp/rc);
1278         Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1279         if(rlim <= rc){
1280           accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1281           accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1282         }
1283         if(rlim > rc){
1284           accetamax = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1285           accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1286             ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1287             rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1288         }
1289       }
1290
1291       if(detamin < eMin && dphimax > pMax ){
1292         Float_t he = eMax + etaJet[k];
1293         Float_t hp = pMax - phiJet[k];
1294         Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1295         Float_t alphae = TMath::ACos(he/rc);
1296         Float_t alphap = TMath::ACos(hp/rc);
1297         Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1298         if(rlim <= rc){
1299           accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1300           accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1301         }
1302         if(rlim > rc){
1303           accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1304           accphimax = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1305             ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1306             rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1307         }
1308       }
1309
1310       if(detamin < eMin && dphimin < pMin ){
1311         Float_t he = eMax + etaJet[k];
1312         Float_t hp = phiJet[k] - pMin;
1313         Float_t rlim = TMath::Sqrt(pow(he,2)+pow(hp,2));
1314         Float_t alphae = TMath::ACos(he/rc);
1315         Float_t alphap = TMath::ACos(hp/rc);
1316         Float_t alphad = (alphae+alphap)/2-TMath::Pi()/4;
1317         if(rlim <= rc){
1318           accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1319           accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp);
1320         }
1321         if(rlim > rc){
1322           accetamin = rc*rc*alphae - he*TMath::Sqrt(rc*rc - he*he);
1323           accphimin = rc*rc*alphap - hp*TMath::Sqrt(rc*rc - hp*hp)-
1324             ((TMath::Sqrt(pow(rc,2)-pow(he,2))-hp)*(TMath::Sqrt(pow(rc,2)-pow(hp,2))-he))/2+
1325             rc*rc*alphad - rc*rc*TMath::Sin(alphad)*TMath::Cos(alphad);
1326         }
1327       }
1328       areaJet[k] = rc*rc*TMath::Pi() - accetamax - accetamin - accphimax - accphimin;
1329       areaOut = areaOut - areaJet[k];
1330     } // end loop on jets
1331
1332     //subtract background using area method
1333     for(Int_t ljet=0; ljet<nJ; ljet++){
1334       Float_t areaRatio = areaJet[ljet]/areaOut;
1335       etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1336     }
1337
1338     // estimate new total background
1339     Float_t areaT = 2*(header->GetLegoEtaMax()*header->GetLegoPhiMax());
1340     etbgTotalN = etOut*areaT/areaOut;
1341   }
1342
1343 }
1344
1345 ////////////////////////////////////////////////////////////////////////
1346 void AliUA1JetFinderV2::SubtractBackgC(const Int_t& nIn, const Int_t&nJ, Float_t&etbgTotalN,
1347                                        const Float_t* ptT, const Float_t* etaT, const Float_t* phiT,
1348                                        Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1349                                        Int_t* multJet, Int_t* injet)
1350 {
1351   //background subtraction using cone method but without correction in dE/deta distribution
1352
1353   //calculate energy inside and outside cones
1354   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1355   Float_t rc= header->GetRadius();
1356   Float_t etIn[30];
1357   Float_t etOut = 0;
1358   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1359     // if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
1360     for(Int_t ijet=0; ijet<nJ; ijet++){
1361       Float_t deta = etaT[jpart] - etaJet[ijet];
1362       Float_t dphi = phiT[jpart] - phiJet[ijet];
1363       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1364       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1365       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1366       if(dr <= rc){ // particles inside this cone
1367         multJet[ijet]++;
1368         injet[jpart] = ijet;
1369         if((fReader->GetCutFlag(jpart)) == 1){ // pt cut
1370           etIn[ijet] += ptT[jpart];
1371           if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
1372         }
1373         break;
1374       }
1375     }// end jets loop
1376     if(injet[jpart] == -1 && fReader->GetCutFlag(jpart) == 1)
1377       etOut += ptT[jpart]; // particle outside cones and pt cut
1378   } //end particle loop
1379
1380   //estimate jets and background areas
1381   Float_t areaJet[30];
1382   Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1383   for(Int_t k=0; k<nJ; k++){
1384     Float_t detamax = etaJet[k] + rc;
1385     Float_t detamin = etaJet[k] - rc;
1386     Float_t accmax = 0.0; Float_t accmin = 0.0;
1387     if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1388       Float_t h = header->GetLegoEtaMax() - etaJet[k];
1389       accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1390     }
1391     if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1392       Float_t h = header->GetLegoEtaMax() + etaJet[k];
1393       accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1394     }
1395     areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1396     areaOut = areaOut - areaJet[k];
1397   }
1398   //subtract background using area method
1399   for(Int_t ljet=0; ljet<nJ; ljet++){
1400     Float_t areaRatio = areaJet[ljet]/areaOut;
1401     etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
1402   }
1403
1404   // estimate new total background
1405   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1406   etbgTotalN = etOut*areaT/areaOut;
1407
1408 }
1409
1410
1411 ////////////////////////////////////////////////////////////////////////
1412 void AliUA1JetFinderV2::SubtractBackgStat(const Int_t& nIn, const Int_t&nJ,Float_t&etbgTotalN,
1413                       const Float_t* ptT, const Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT,
1414                       const Float_t* sFlagT, Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet,
1415                       Float_t* etsigJet, Int_t* multJet, Int_t* injet)
1416 {
1417
1418   //background subtraction using statistical method
1419   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1420   Float_t etbgStat = header->GetBackgStat(); // pre-calculated background
1421
1422   //calculate energy inside
1423   Float_t rc= header->GetRadius();
1424   Float_t etIn[30];
1425
1426   for(Int_t jpart = 0; jpart < nIn; jpart++)
1427     { // loop for all particles in array
1428
1429       for(Int_t ijet=0; ijet<nJ; ijet++)
1430         {
1431           Float_t deta = etaT[jpart] - etaJet[ijet];
1432           Float_t dphi = phiT[jpart] - phiJet[ijet];
1433           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1434           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1435           Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1436           if(dr <= rc){ // particles inside this cone
1437             multJet[ijet]++;
1438             injet[jpart] = ijet;
1439             if(cFlagT[jpart] == 1){ // pt cut
1440               etIn[ijet]+= ptT[jpart];
1441               if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1442             }
1443             break;
1444           }
1445         }// end jets loop
1446     } //end particle loop
1447
1448   //calc jets areas
1449   Float_t areaJet[30];
1450   Float_t areaOut = 4*(header->GetLegoEtaMax())*TMath::Pi();
1451   for(Int_t k=0; k<nJ; k++)
1452     {
1453       Float_t detamax = etaJet[k] + rc;
1454       Float_t detamin = etaJet[k] - rc;
1455       Float_t accmax = 0.0; Float_t accmin = 0.0;
1456       if(detamax > header->GetLegoEtaMax()){ // sector outside etamax
1457         Float_t h = header->GetLegoEtaMax() - etaJet[k];
1458         accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1459       }
1460       if(detamin < header->GetLegoEtaMin()){ // sector outside etamin
1461       Float_t h = header->GetLegoEtaMax() + etaJet[k];
1462       accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
1463       }
1464       areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
1465     }
1466
1467   //subtract background using area method
1468   for(Int_t ljet=0; ljet<nJ; ljet++){
1469     Float_t areaRatio = areaJet[ljet]/areaOut;
1470     etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
1471   }
1472
1473   etbgTotalN = etbgStat;
1474 }
1475
1476 ////////////////////////////////////////////////////////////////////////
1477 void AliUA1JetFinderV2::SubtractBackgCone(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN, Float_t* ptT,
1478                                           Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
1479                                           Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1480                                           Int_t* multJet, Int_t* injet)
1481 {
1482   // Cone background subtraction method taking into acount dEt/deta distribution
1483   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1484   //general
1485   Float_t rc= header->GetRadius();
1486   Float_t etamax = header->GetLegoEtaMax();
1487   Float_t etamin = header->GetLegoEtaMin();
1488   Int_t ndiv = 100;
1489
1490   // jet energy and area arrays
1491   TH1F* hEtJet[30];
1492   TH1F* hAreaJet[30];
1493   for(Int_t mjet=0; mjet<nJ; mjet++){
1494     char hEtname[256]; char hAreaname[256];
1495     sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
1496     hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
1497     hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
1498   }
1499   // background energy and area
1500   TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
1501   TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
1502
1503   //fill energies
1504   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1505     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
1506       Float_t deta = etaT[jpart] - etaJet[ijet];
1507       Float_t dphi = phiT[jpart] - phiJet[ijet];
1508       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1509       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1510       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1511       if(dr <= rc){ // particles inside this cone
1512         injet[jpart] = ijet;
1513         multJet[ijet]++;
1514         if(cFlagT[jpart] == 1){// pt cut
1515           hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
1516           if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1517         }
1518         break;
1519       }
1520     }// end jets loop
1521
1522     if(injet[jpart] == -1  && cFlagT[jpart] == 1)
1523       hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
1524   } //end particle loop
1525
1526   //calc areas
1527   Float_t eta0 = etamin;
1528   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
1529   Float_t eta1 = eta0 + etaw;
1530   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
1531     Float_t etac = eta0 + etaw/2.0;
1532     Float_t areabg = etaw*2.0*TMath::Pi();
1533     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
1534       Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
1535       Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
1536       Float_t acc0 = 0.0; Float_t acc1 = 0.0;
1537       Float_t areaj = 0.0;
1538       if(deta0 > rc && deta1 < rc){
1539         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1540         areaj = acc1;
1541       }
1542       if(deta0 < rc && deta1 > rc){
1543         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1544         areaj = acc0;
1545       }
1546       if(deta0 < rc && deta1 < rc){
1547         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1548         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1549         if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
1550         if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
1551         if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
1552       }
1553       hAreaJet[ijet]->Fill(etac,areaj);
1554       areabg = areabg - areaj;
1555     } // end jets loop
1556     hAreaBackg->Fill(etac,areabg);
1557     eta0 = eta1;
1558     eta1 = eta1 + etaw;
1559   } // end loop for all eta bins
1560
1561   //subtract background
1562   for(Int_t kjet=0; kjet<nJ; kjet++){
1563     etJet[kjet] = 0.0; // first  clear etJet for this jet
1564     for(Int_t bin = 0; bin< ndiv; bin++){
1565       if(hAreaJet[kjet]->GetBinContent(bin)){
1566         Float_t areab = hAreaBackg->GetBinContent(bin);
1567         Float_t etb = hEtBackg->GetBinContent(bin);
1568         Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
1569         etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
1570       }
1571     }
1572   }
1573
1574   // calc background total
1575   Double_t etOut = hEtBackg->Integral();
1576   Double_t areaOut = hAreaBackg->Integral();
1577   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1578   etbgTotalN = etOut*areaT/areaOut;
1579
1580   //delete
1581   for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
1582     delete hEtJet[ljet];
1583     delete hAreaJet[ljet];
1584   }
1585
1586   delete hEtBackg;
1587   delete hAreaBackg;
1588 }
1589
1590 ////////////////////////////////////////////////////////////////////////
1591 void AliUA1JetFinderV2::SubtractBackgRatio(const Int_t& nIn, const Int_t&nJ,Float_t& etbgTotalN,
1592                       Float_t* ptT, Float_t* etaT, const Float_t* phiT, const Float_t* cFlagT, const Float_t* sFlagT,
1593                       Float_t* etJet, const Float_t* etaJet, const Float_t* phiJet, Float_t* etsigJet,
1594                       Int_t* multJet, Int_t* injet)
1595 {
1596   // Ratio background subtraction method taking into acount dEt/deta distribution
1597   AliUA1JetHeaderV1* header = (AliUA1JetHeaderV1*) fHeader;
1598   //factor F calc before
1599   Float_t bgRatioCut = header->GetBackgCutRatio();
1600
1601   //general
1602   Float_t rc= header->GetRadius();
1603   Float_t etamax = header->GetLegoEtaMax();
1604   Float_t etamin = header->GetLegoEtaMin();
1605   Int_t ndiv = 100;
1606
1607   // jet energy and area arrays
1608   TH1F* hEtJet[30];
1609   TH1F* hAreaJet[30];
1610   for(Int_t mjet=0; mjet<nJ; mjet++){
1611     char hEtname[256]; char hAreaname[256];
1612     sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
1613     hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);        // change range
1614     hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);  // change range
1615   }
1616   // background energy and area
1617   TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);         // change range
1618   TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);  // change range
1619
1620   //fill energies
1621   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
1622     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
1623       Float_t deta = etaT[jpart] - etaJet[ijet];
1624       Float_t dphi = phiT[jpart] - phiJet[ijet];
1625       if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
1626       if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
1627       Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
1628       if(dr <= rc){ // particles inside this cone
1629         multJet[ijet]++;
1630         injet[jpart] = ijet;
1631         if(cFlagT[jpart] == 1){ //pt cut
1632           hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
1633           if(sFlagT[jpart] == 1) etsigJet[ijet] += ptT[jpart];
1634         }
1635         break;
1636       }
1637     }// end jets loop
1638     if(injet[jpart] == -1) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
1639   } //end particle loop
1640
1641   //calc areas
1642   Float_t eta0 = etamin;
1643   Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
1644   Float_t eta1 = eta0 + etaw;
1645   for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
1646     Float_t etac = eta0 + etaw/2.0;
1647     Float_t areabg = etaw*2.0*TMath::Pi();
1648     for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
1649       Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
1650       Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
1651       Float_t acc0 = 0.0; Float_t acc1 = 0.0;
1652       Float_t areaj = 0.0;
1653       if(deta0 > rc && deta1 < rc){
1654         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1655         areaj = acc1;
1656       }
1657       if(deta0 < rc && deta1 > rc){
1658         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1659         areaj = acc0;
1660       }
1661       if(deta0 < rc && deta1 < rc){
1662         acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
1663         acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
1664         if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
1665         if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
1666         if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
1667       }
1668       hAreaJet[ijet]->Fill(etac,areaj);
1669       areabg = areabg - areaj;
1670     } // end jets loop
1671     hAreaBackg->Fill(etac,areabg);
1672     eta0 = eta1;
1673     eta1 = eta1 + etaw;
1674   } // end loop for all eta bins
1675
1676   //subtract background
1677   for(Int_t kjet=0; kjet<nJ; kjet++){
1678     etJet[kjet] = 0.0; // first  clear etJet for this jet
1679     for(Int_t bin = 0; bin< ndiv; bin++){
1680       if(hAreaJet[kjet]->GetBinContent(bin)){
1681         Float_t areab = hAreaBackg->GetBinContent(bin);
1682         Float_t etb = hEtBackg->GetBinContent(bin);
1683         Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
1684         etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
1685       }
1686     }
1687   }
1688
1689   // calc background total
1690   Double_t etOut = hEtBackg->Integral();
1691   Double_t areaOut = hAreaBackg->Integral();
1692   Float_t areaT = 4*(header->GetLegoEtaMax())*TMath::Pi();
1693   etbgTotalN = etOut*areaT/areaOut;
1694
1695   //delete
1696   for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
1697     delete hEtJet[ljet];
1698     delete hAreaJet[ljet];
1699   }
1700
1701   delete hEtBackg;
1702   delete hAreaBackg;
1703 }
1704
1705 ////////////////////////////////////////////////////////////////////////
1706 void AliUA1JetFinderV2::Reset()
1707 {
1708   fLego->Reset();
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 }