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