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