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