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