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