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