New version of the jet finder (Rafael)
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinderV1.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 //---------------------------------------------------------------------
17 // UA1 Cone Algorithm Jet finder
18 // manages the search for jets
19 // Author: Rafael.Diaz.Valdes@cern.ch
20 // (version in c++)
21 //---------------------------------------------------------------------
22
23 #include <TLorentzVector.h>
24 #include <TFile.h>
25 #include <TH1F.h>
26 #include <TH2F.h>
27 #include <TArrayF.h>
28 #include "AliUA1JetFinderV1.h"
29 #include "AliUA1JetHeaderV1.h"
30 #include "AliJetReaderHeader.h"
31 #include "AliJetReader.h"
32 #include "AliJet.h"
33
34
35 ClassImp(AliUA1JetFinderV1)
36
37 ////////////////////////////////////////////////////////////////////////
38
39 AliUA1JetFinderV1::AliUA1JetFinderV1()
40
41 {
42   // Constructor
43   fHeader = 0x0;
44   fLego   = 0x0;
45 }
46
47 ////////////////////////////////////////////////////////////////////////
48
49 AliUA1JetFinderV1::~AliUA1JetFinderV1()
50
51 {
52   // destructor
53 }
54
55 ////////////////////////////////////////////////////////////////////////
56
57
58 void AliUA1JetFinderV1::FindJets()
59
60 {
61   //1) Fill cell map array
62   //2) calculate total energy and fluctuation level
63   //3) Run algorithm
64   //   3.1) look centroides in cell map
65   //   3.2) calculate total energy in cones
66   //   3.3) flag as a possible jet
67   //   3.4) reorder cones by energy
68   //4) subtract backg in accepted jets
69   //5) fill AliJet list
70
71   // transform input to pt,eta,phi plus lego
72   TClonesArray *lvArray = fReader->GetMomentumArray();
73   Int_t nIn =  lvArray->GetEntries();
74   if (nIn == 0) return;
75
76   // local arrays for input
77   Float_t* ptT  = new Float_t[nIn];
78   Float_t* etaT = new Float_t[nIn];
79   Float_t* phiT = new Float_t[nIn];
80   Int_t*   injet = new Int_t[nIn];
81
82   //total energy in array
83   Float_t  etbgTotal = 0.0;
84   TH1F* hPtTotal = new TH1F("hPt","Pt distribution of all particles ",100,0.0,15.0);
85
86   // load input vectors and calculate total energy in array
87   for (Int_t i = 0; i < nIn; i++){
88     TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
89     ptT[i]  = lv->Pt();
90     etaT[i] = lv->Eta();
91     phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
92     if (fReader->GetCutFlag(i) != 1) continue;
93     fLego->Fill(etaT[i], phiT[i], ptT[i]);
94     hPtTotal->Fill(ptT[i]);
95     etbgTotal+= ptT[i];
96   }
97   fJets->SetNinput(nIn);
98
99   // calculate total energy and fluctuation in map
100   Double_t meanpt = hPtTotal->GetMean();
101   Double_t ptRMS = hPtTotal->GetRMS();
102   Double_t npart = hPtTotal->GetEntries();
103   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
104
105   // arrays to hold jets
106   Float_t* etaJet = new Float_t[30];
107   Float_t* phiJet = new Float_t[30];
108   Float_t* etJet  = new Float_t[30];
109   Float_t* etsigJet  = new Float_t[30]; //signal et in jet
110   Float_t* etallJet = new Float_t[30];  // total et in jet (tmp variable)
111   Int_t* ncellsJet = new Int_t[30];
112   Int_t* multJet  = new Int_t[30];
113   Int_t nJets; // to hold number of jets found by algorithm
114   Int_t nj;    // number of jets accepted
115   Float_t prec  = fHeader->GetPrecBg();
116   Float_t bgprec = 1;
117   while(bgprec > prec){
118      //reset jet arrays in memory
119      memset(etaJet,0,sizeof(Float_t)*30);
120      memset(phiJet,0,sizeof(Float_t)*30);
121      memset(etJet,0,sizeof(Float_t)*30);
122      memset(etallJet,0,sizeof(Float_t)*30);
123      memset(etsigJet,0,sizeof(Float_t)*30);
124      memset(ncellsJet,0,sizeof(Int_t)*30);
125      memset(multJet,0,sizeof(Int_t)*30);
126      nJets = 0;
127      nj = 0;
128      // reset particles-jet array in memory
129      memset(injet,0,sizeof(Int_t)*nIn);
130      //run cone algorithm finder
131      RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
132      //run background subtraction
133      if(nJets > fHeader->GetNAcceptJets()) // limited number of accepted jets per event
134        nj = fHeader->GetNAcceptJets();
135      else
136        nj = nJets;
137      //subtract background
138      Float_t etbgTotalN = 0.0; //new background
139      if(fHeader->GetBackgMode() == 1) // standar
140         SubtractBackg(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
141      if(fHeader->GetBackgMode() == 2) //cone
142         SubtractBackgCone(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
143      if(fHeader->GetBackgMode() == 3) //ratio
144         SubtractBackgRatio(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
145      if(fHeader->GetBackgMode() == 4) //statistic
146         SubtractBackgStat(nIn,nj,etbgTotalN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJet,injet);
147      //calc precision
148      if(etbgTotalN != 0.0)
149         bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
150      else
151         bgprec = 0;
152      etbgTotal = etbgTotalN; // update with new background estimation
153   } //end while
154
155   // add jets to list
156   Int_t* idxjets = new Int_t[nj];
157   Int_t nselectj = 0;
158   for(Int_t kj=0; kj<nj; kj++){
159      if(etJet[kj] > fHeader->GetMinJetEt()){ // check if et jets > et min
160       Float_t px, py,pz,en; // convert to 4-vector
161       px = etJet[kj] * TMath::Cos(phiJet[kj]);
162       py = etJet[kj] * TMath::Sin(phiJet[kj]);
163       pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj])));
164       en = TMath::Sqrt(px * px + py * py + pz * pz);
165       fJets->AddJet(px, py, pz, en);
166       idxjets[nselectj] = kj;
167       nselectj++;
168      }
169   }
170   //add signal percentage and total signal  in AliJets for analysis tool
171   Float_t* percentage  = new Float_t[nselectj];
172   Int_t* ncells      = new Int_t[nselectj];
173   Int_t* mult        = new Int_t[nselectj];
174
175   for(Int_t i = 0; i< nselectj; i++){
176      percentage[i] = etsigJet[idxjets[i]]/etJet[idxjets[i]];
177      ncells[i] = ncellsJet[idxjets[i]];
178      mult[i] = multJet[idxjets[i]];
179   }
180   fJets->SetNCells(ncells);
181   fJets->SetPtFromSignal(percentage);
182   fJets->SetMultiplicities(mult);
183   fJets->SetInJet(injet);
184   fJets->SetEtaIn(etaT);
185   fJets->SetPhiIn(phiT);
186   fJets->SetPtIn(ptT);
187   fJets->SetEtAvg(etbgTotal/(4*(fHeader->GetLegoEtaMax())*TMath::Pi()));
188
189
190   //delete
191   delete ptT;
192   delete etaT;
193   delete phiT;
194   delete injet;
195   delete hPtTotal;
196   delete etaJet;
197   delete phiJet;
198   delete etJet;
199   delete etsigJet;
200   delete etallJet;
201   delete ncellsJet;
202   delete multJet;
203   delete idxjets;
204   delete percentage;
205   delete ncells;
206   delete mult;
207
208
209 }
210
211 ////////////////////////////////////////////////////////////////////////
212
213 void AliUA1JetFinderV1::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
214                                   Float_t* etJet,Float_t* etaJet, Float_t* phiJet,
215                                   Float_t* etallJet, Int_t* ncellsJet)
216 {
217
218    //dump lego
219   // check enough space! *to be done*
220   Float_t etCell[60000];   //! Cell Energy
221   Float_t etaCell[60000];  //! Cell eta
222   Float_t phiCell[60000];  //! Cell phi
223   Int_t   flagCell[60000]; //! Cell flag
224
225   Int_t nCell = 0;
226   TAxis* xaxis = fLego->GetXaxis();
227   TAxis* yaxis = fLego->GetYaxis();
228   Float_t e = 0.0;
229   for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
230       for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
231                e = fLego->GetBinContent(i,j);
232                if (e < 0.0) continue; // don't include this cells
233                Float_t eta  = xaxis->GetBinCenter(i);
234                Float_t phi  = yaxis->GetBinCenter(j);
235                etCell[nCell]  = e;
236                etaCell[nCell] = eta;
237                phiCell[nCell] = phi;
238           flagCell[nCell] = 0; //default
239                nCell++;
240       }
241   }
242
243   // Parameters from header
244   Float_t minmove = fHeader->GetMinMove();
245   Float_t maxmove = fHeader->GetMaxMove();
246   Float_t rc      = fHeader->GetRadius();
247   Float_t etseed  = fHeader->GetEtSeed();
248   //Float_t etmin   = fHeader->GetMinJetEt();
249
250
251
252   // tmp array of jets form algoritm
253   Float_t etaAlgoJet[30];
254   Float_t phiAlgoJet[30];
255   Float_t etAlgoJet[30];
256   Int_t   ncellsAlgoJet[30];
257
258   //run algorithm//
259
260   // sort cells by et
261   Int_t * index  = new Int_t[nCell];
262   TMath::Sort(nCell, etCell, index);
263   // variable used in centroide loop
264   Float_t eta = 0.0;
265   Float_t phi = 0.0;
266   Float_t eta0 = 0.0;
267   Float_t phi0 = 0.0;
268   Float_t etab = 0.0;
269   Float_t phib = 0.0;
270   Float_t etas = 0.0;
271   Float_t phis = 0.0;
272   Float_t ets = 0.0;
273   Float_t deta = 0.0;
274   Float_t dphi = 0.0;
275   Float_t dr = 0.0;
276   Float_t etsb = 0.0;
277   Float_t etasb = 0.0;
278   Float_t phisb = 0.0;
279
280
281   for(Int_t icell = 0; icell < nCell; icell++){
282         Int_t jcell = index[icell];
283         if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
284         if(flagCell[jcell] != 0) continue; // if cell was used before
285         eta  = etaCell[jcell];
286         phi  = phiCell[jcell];
287         eta0 = eta;
288         phi0 = phi;
289         etab = eta;
290         phib = phi;
291         ets  = etCell[jcell];
292         etas = 0.0;
293         phis = 0.0;
294         etsb = ets;
295         etasb = 0.0;
296         phisb = 0.0;
297         for(Int_t kcell =0; kcell < nCell; kcell++){
298             Int_t lcell = index[kcell];
299             if(lcell == jcell) continue; // cell itself
300             if(flagCell[lcell] != 0) continue; // cell used before
301             if(etCell[lcell] > etCell[jcell]) continue;
302             //calculate dr
303             deta = etaCell[lcell] - eta;
304                  dphi = phiCell[lcell] - phi;
305                  if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
306                  if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
307                  dr = TMath::Sqrt(deta * deta + dphi * dphi);
308             if(dr <= rc){
309                // calculate offset from initiate cell
310                deta = etaCell[lcell] - eta0;
311                dphi = phiCell[lcell] - phi0;
312                if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
313                     if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
314                etas = etas + etCell[lcell]*deta;
315                phis = phis + etCell[lcell]*dphi;
316                ets = ets + etCell[lcell];
317                //new weighted eta and phi including this cell
318                eta = eta0 + etas/ets;
319                phi = phi0 + phis/ets;
320                // if cone does not move much, just go to next step
321                dr = TMath::Sqrt((eta-etab)*(eta-etab) + (phi-phib)*(phi-phib));
322                if(dr <= minmove) break;
323                // cone should not move more than max_mov
324                dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
325                if(dr > maxmove){
326                   eta = etab;
327                   phi = phib;
328                   ets = etsb;
329                   etas = etasb;
330                   phis = phisb;
331                }else{ // store this loop information
332                  etab=eta;
333                  phib=phi;
334                  etsb = ets;
335                  etasb = etas;
336                  phisb = phis;
337                }
338             }
339         }//end of cells loop looking centroide
340
341         //avoid cones overloap (to be implemented in the future)
342
343         //flag cells in Rc, estimate total energy in cone
344         Float_t etCone = 0.0;
345         Int_t   nCellIn = 0;
346         rc = fHeader->GetRadius();
347         for(Int_t ncell =0; ncell < nCell; ncell++){
348             if(flagCell[ncell] != 0) continue; // cell used before
349            //calculate dr
350             deta = etaCell[ncell] - eta;
351                  dphi = phiCell[ncell] - phi;
352                  if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
353                  if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
354                  dr = TMath::Sqrt(deta * deta + dphi * dphi);
355             if(dr <= rc){  // cell in cone
356                flagCell[ncell] = -1;
357                etCone+=etCell[ncell];
358                nCellIn++;
359             }
360         }
361
362         // select jets with et > background
363         // estimate max fluctuation of background in cone
364         Double_t ncellin = (Double_t)nCellIn;
365         Double_t ntcell  = (Double_t)nCell;
366         Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
367         // min cone et
368         Double_t etcmin = etCone ;  // could be used etCone - etmin !!
369         //desicions !! etbmax < etcmin
370         for(Int_t mcell =0; mcell < nCell; mcell++){
371             if(flagCell[mcell] == -1){
372               if(etbmax < etcmin)
373                  flagCell[mcell] = 1; //flag cell as used
374               else
375                  flagCell[mcell] = 0; // leave it free
376             }
377         }
378         //store tmp jet info !!!
379         // reject tmp outside acceptable eta range
380        if ((eta < (fHeader->GetJetEtaMax())) &&
381            (eta > (fHeader->GetJetEtaMin())) &&
382            (etbmax < etcmin) ){
383              etaAlgoJet[nJets] = eta;
384              phiAlgoJet[nJets] = phi;
385              etAlgoJet[nJets] = etCone;
386              ncellsAlgoJet[nJets] = nCellIn;
387              nJets++;
388         }
389
390   } // end of cells loop
391
392   //reorder jets by et in cone
393   //sort jets by energy
394   Int_t * idx  = new Int_t[nJets];
395   TMath::Sort(nJets, etAlgoJet, idx);
396   for(Int_t p = 0; p < nJets; p++){
397      etaJet[p] = etaAlgoJet[idx[p]];
398      phiJet[p] = phiAlgoJet[idx[p]];
399      etJet[p] = etAlgoJet[idx[p]];
400      etallJet[p] = etAlgoJet[idx[p]];
401      ncellsJet[p] = ncellsAlgoJet[idx[p]];
402   }
403
404
405   //delete
406   delete index;
407   delete idx;
408
409 }
410 ////////////////////////////////////////////////////////////////////////
411
412 void AliUA1JetFinderV1::SubtractBackg(Int_t& nIn, Int_t&nJ, Float_t&etbgTotalN,
413                       Float_t* ptT, Float_t* etaT, Float_t* phiT,
414                       Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
415                       Int_t* multJet, Int_t* injet)
416 {
417   //background subtraction using cone method but without correction in dE/deta distribution
418
419   //calculate energy inside and outside cones
420   Float_t rc= fHeader->GetRadius();
421   Float_t etIn[30];
422   Float_t etOut = 0;
423   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
424      if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
425      for(Int_t ijet=0; ijet<nJ; ijet++){
426          Float_t deta = etaT[jpart] - etaJet[ijet];
427               Float_t dphi = phiT[jpart] - phiJet[ijet];
428          if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
429               if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
430               Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
431          if(dr <= rc){ // particles inside this cone
432              etIn[ijet] += ptT[jpart];
433              if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet]+= ptT[jpart];
434              multJet[ijet]++;
435              injet[jpart] = ijet;
436              break;
437          }
438      }// end jets loop
439      if(injet[jpart] == 0) etOut += ptT[jpart]; // particle outside cones
440   } //end particle loop
441
442   //estimate jets and background areas
443   Float_t areaJet[30];
444   Float_t areaOut = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
445   for(Int_t k=0; k<nJ; k++){
446       Float_t detamax = etaJet[k] + rc;
447       Float_t detamin = etaJet[k] - rc;
448       Float_t accmax = 0.0; Float_t accmin = 0.0;
449       if(detamax > fHeader->GetLegoEtaMax()){ // sector outside etamax
450          Float_t h = fHeader->GetLegoEtaMax() - etaJet[k];
451          accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
452       }
453       if(detamin < fHeader->GetLegoEtaMin()){ // sector outside etamin
454          Float_t h = fHeader->GetLegoEtaMax() + etaJet[k];
455          accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
456       }
457       areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
458       areaOut = areaOut - areaJet[k];
459   }
460   //subtract background using area method
461   for(Int_t ljet=0; ljet<nJ; ljet++){
462      Float_t areaRatio = areaJet[ljet]/areaOut;
463      etJet[ljet] = etIn[ljet]-etOut*areaRatio; // subtraction
464   }
465
466   // estimate new total background
467   Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
468   etbgTotalN = etOut*areaT/areaOut;
469
470
471 }
472
473 ////////////////////////////////////////////////////////////////////////
474
475 void AliUA1JetFinderV1::SubtractBackgStat(Int_t& nIn, Int_t&nJ,Float_t&etbgTotalN,
476                       Float_t* ptT, Float_t* etaT, Float_t* phiT,
477                       Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
478                       Int_t* multJet, Int_t* injet)
479 {
480
481   //background subtraction using statistical method
482
483   Float_t etbgStat = fHeader->GetBackgStat(); // pre-calculated background
484
485   //calculate energy inside
486   Float_t rc= fHeader->GetRadius();
487   Float_t etIn[30];
488
489   for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
490      if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
491      for(Int_t ijet=0; ijet<nJ; ijet++){
492          Float_t deta = etaT[jpart] - etaJet[ijet];
493               Float_t dphi = phiT[jpart] - phiJet[ijet];
494          if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
495               if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
496               Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
497          if(dr <= rc){ // particles inside this cone
498              etIn[ijet]+= ptT[jpart];
499              if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
500              multJet[ijet]++;
501              injet[jpart] = ijet;
502              break;
503          }
504      }// end jets loop
505   } //end particle loop
506
507   //calc jets areas
508   Float_t areaJet[30];
509   Float_t areaOut = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
510   for(Int_t k=0; k<nJ; k++){
511       Float_t detamax = etaJet[k] + rc;
512       Float_t detamin = etaJet[k] - rc;
513       Float_t accmax = 0.0; Float_t accmin = 0.0;
514       if(detamax > fHeader->GetLegoEtaMax()){ // sector outside etamax
515          Float_t h = fHeader->GetLegoEtaMax() - etaJet[k];
516          accmax = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
517       }
518       if(detamin < fHeader->GetLegoEtaMin()){ // sector outside etamin
519          Float_t h = fHeader->GetLegoEtaMax() + etaJet[k];
520          accmin = rc*rc*TMath::ACos(h/rc) - h*TMath::Sqrt(rc*rc - h*h);
521       }
522       areaJet[k] = rc*rc*TMath::Pi() - accmax - accmin;
523   }
524
525   //subtract background using area method
526   for(Int_t ljet=0; ljet<nJ; ljet++){
527      Float_t areaRatio = areaJet[ljet]/areaOut;
528      etJet[ljet] = etIn[ljet]-etbgStat*areaRatio; // subtraction
529   }
530
531   etbgTotalN = etbgStat;
532
533 }
534
535 ////////////////////////////////////////////////////////////////////////
536
537 void AliUA1JetFinderV1::SubtractBackgCone(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN,
538                       Float_t* ptT, Float_t* etaT, Float_t* phiT,
539                       Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
540                       Int_t* multJet, Int_t* injet)
541 {
542    // Cone background subtraction method taking into acount dEt/deta distribution
543
544    //general
545    Float_t rc= fHeader->GetRadius();
546    Float_t etamax = fHeader->GetLegoEtaMax();
547    Float_t etamin = fHeader->GetLegoEtaMin();
548    Int_t ndiv = 100;
549
550    // jet energy and area arrays
551    TH1F* hEtJet[30];
552    TH1F* hAreaJet[30];
553    for(Int_t mjet=0; mjet<nJ; mjet++){
554      char hEtname[256]; char hAreaname[256];
555      sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
556      hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);
557      hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);
558   }
559    // background energy and area
560    TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);
561    TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);
562
563    //fill energies
564    for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
565      if((fReader->GetCutFlag(jpart)) != 1) continue; // pt cut
566      for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
567          Float_t deta = etaT[jpart] - etaJet[ijet];
568               Float_t dphi = phiT[jpart] - phiJet[ijet];
569          if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
570               if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
571               Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
572          if(dr <= rc){ // particles inside this cone
573              hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone
574              multJet[ijet]++;
575              if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
576              injet[jpart] = ijet;
577              break;
578          }
579      }// end jets loop
580      if(injet[jpart] == 0) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
581   } //end particle loop
582
583    //calc areas
584    Float_t eta0 = etamin;
585    Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
586    Float_t eta1 = eta0 + etaw;
587    for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
588       Float_t etac = eta0 + etaw/2.0;
589       Float_t areabg = etaw*2.0*TMath::Pi();
590       for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
591           Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
592           Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
593           Float_t acc0 = 0.0; Float_t acc1 = 0.0;
594           Float_t areaj = 0.0;
595           if(deta0 > rc && deta1 < rc){
596              acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
597              areaj = acc1;
598           }
599           if(deta0 < rc && deta1 > rc){
600              acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
601              areaj = acc0;
602           }
603           if(deta0 < rc && deta1 < rc){
604              acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
605              acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
606              if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
607              if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
608              if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
609           }
610           hAreaJet[ijet]->Fill(etac,areaj);
611           areabg = areabg - areaj;
612       } // end jets loop
613       hAreaBackg->Fill(etac,areabg);
614       eta0 = eta1;
615       eta1 = eta1 + etaw;
616    } // end loop for all eta bins
617
618    //subtract background
619    for(Int_t kjet=0; kjet<nJ; kjet++){
620        etJet[kjet] = 0.0; // first  clear etJet for this jet
621        for(Int_t bin = 0; bin< ndiv; bin++){
622            if(hAreaJet[kjet]->GetBinContent(bin)){
623               Float_t areab = hAreaBackg->GetBinContent(bin);
624               Float_t etb = hEtBackg->GetBinContent(bin);
625               Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
626               etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR); //subtraction
627            }
628        }
629    }
630
631    // calc background total
632    Double_t etOut = hEtBackg->Integral();
633    Double_t areaOut = hAreaBackg->Integral();
634    Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
635    etbgTotalN = etOut*areaT/areaOut;
636
637    //delete
638    for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
639        delete hEtJet[ljet];
640        delete hAreaJet[ljet];
641    }
642
643    delete hEtBackg;
644    delete hAreaBackg;
645 }
646
647 ////////////////////////////////////////////////////////////////////////
648
649
650 void AliUA1JetFinderV1::SubtractBackgRatio(Int_t& nIn, Int_t&nJ,Float_t& etbgTotalN,
651                       Float_t* ptT, Float_t* etaT, Float_t* phiT,
652                       Float_t* etJet,Float_t* etaJet, Float_t* phiJet, Float_t* etsigJet,
653                        Int_t* multJet, Int_t* injet)
654 {
655    // Ratio background subtraction method taking into acount dEt/deta distribution
656
657    //factor F calc before
658     Float_t bgRatioCut = fHeader->GetBackgCutRatio();
659
660
661    //general
662    Float_t rc= fHeader->GetRadius();
663    Float_t etamax = fHeader->GetLegoEtaMax();
664    Float_t etamin = fHeader->GetLegoEtaMin();
665    Int_t ndiv = 100;
666
667    // jet energy and area arrays
668    TH1F* hEtJet[30];
669    TH1F* hAreaJet[30];
670    for(Int_t mjet=0; mjet<nJ; mjet++){
671      char hEtname[256]; char hAreaname[256];
672      sprintf(hEtname, "hEtJet%d", mjet); sprintf(hAreaname, "hAreaJet%d", mjet);
673      hEtJet[mjet] = new TH1F(hEtname,"et dist in eta ",ndiv,etamin,etamax);        // change range
674      hAreaJet[mjet] = new TH1F(hAreaname,"area dist in eta ",ndiv,etamin,etamax);  // change range
675   }
676    // background energy and area
677    TH1F* hEtBackg = new TH1F("hEtBackg"," backg et dist in eta ",ndiv,etamin,etamax);         // change range
678    TH1F* hAreaBackg = new TH1F("hAreaBackg","backg area dist in eta ",ndiv,etamin,etamax);  // change range
679
680    //fill energies
681    for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
682      if((fReader->GetCutFlag(jpart)) != 1) continue;
683      for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
684          Float_t deta = etaT[jpart] - etaJet[ijet];
685               Float_t dphi = phiT[jpart] - phiJet[ijet];
686          if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
687               if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
688               Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
689          if(dr <= rc){ // particles inside this cone
690             hEtJet[ijet]->Fill(etaT[jpart],ptT[jpart]); //particle inside cone and pt cut
691             if(fReader->GetSignalFlag(jpart) == 1) etsigJet[ijet] += ptT[jpart];
692             multJet[ijet]++;
693             injet[jpart] = ijet;
694             break;
695          }
696      }// end jets loop
697      if(injet[jpart] == 0) hEtBackg->Fill(etaT[jpart],ptT[jpart]); // particle outside cones
698   } //end particle loop
699
700    //calc areas
701    Float_t eta0 = etamin;
702    Float_t etaw = (etamax - etamin)/((Float_t)ndiv);
703    Float_t eta1 = eta0 + etaw;
704    for(Int_t etabin = 0; etabin< ndiv; etabin++){ // loop for all eta bins
705       Float_t etac = eta0 + etaw/2.0;
706       Float_t areabg = etaw*2.0*TMath::Pi();
707       for(Int_t ijet=0; ijet<nJ; ijet++){  // loop for all jets
708           Float_t deta0 = TMath::Abs(eta0 - etaJet[ijet]);
709           Float_t deta1 = TMath::Abs(eta1 - etaJet[ijet]);
710           Float_t acc0 = 0.0; Float_t acc1 = 0.0;
711           Float_t areaj = 0.0;
712           if(deta0 > rc && deta1 < rc){
713              acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
714              areaj = acc1;
715           }
716           if(deta0 < rc && deta1 > rc){
717              acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
718              areaj = acc0;
719           }
720           if(deta0 < rc && deta1 < rc){
721              acc0 = rc*rc*TMath::ACos(deta0/rc) - deta0*TMath::Sqrt(rc*rc - deta0*deta0);
722              acc1 = rc*rc*TMath::ACos(deta1/rc) - deta1*TMath::Sqrt(rc*rc - deta1*deta1);
723              if(eta1<etaJet[ijet]) areaj = acc1-acc0;  // case 1
724              if((eta0 < etaJet[ijet]) && (etaJet[ijet]<eta1)) areaj = rc*rc*TMath::Pi() - acc1 -acc0; // case 2
725              if(etaJet[ijet] < eta0) areaj = acc0 -acc1; // case 3
726           }
727           hAreaJet[ijet]->Fill(etac,areaj);
728           areabg = areabg - areaj;
729       } // end jets loop
730       hAreaBackg->Fill(etac,areabg);
731       eta0 = eta1;
732       eta1 = eta1 + etaw;
733    } // end loop for all eta bins
734
735    //subtract background
736    for(Int_t kjet=0; kjet<nJ; kjet++){
737        etJet[kjet] = 0.0; // first  clear etJet for this jet
738        for(Int_t bin = 0; bin< ndiv; bin++){
739            if(hAreaJet[kjet]->GetBinContent(bin)){
740               Float_t areab = hAreaBackg->GetBinContent(bin);
741               Float_t etb = hEtBackg->GetBinContent(bin);
742               Float_t areaR = (hAreaJet[kjet]->GetBinContent(bin))/areab;
743               etJet[kjet] = etJet[kjet] + ((hEtJet[kjet]->GetBinContent(bin)) - etb*areaR*bgRatioCut); //subtraction
744            }
745        }
746    }
747
748    // calc background total
749    Double_t etOut = hEtBackg->Integral();
750    Double_t areaOut = hAreaBackg->Integral();
751    Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
752    etbgTotalN = etOut*areaT/areaOut;
753
754    //delete
755    for(Int_t ljet=0; ljet<nJ; ljet++){  // loop for all jets
756        delete hEtJet[ljet];
757        delete hAreaJet[ljet];
758    }
759
760    delete hEtBackg;
761    delete hAreaBackg;
762 }
763
764 ////////////////////////////////////////////////////////////////////////
765
766
767 void AliUA1JetFinderV1::Reset()
768 {
769   fLego->Reset();
770   fJets->ClearJets();
771 }
772
773 ////////////////////////////////////////////////////////////////////////
774
775 void AliUA1JetFinderV1::WriteJHeaderToFile()
776 {
777   fOut->cd();
778   fHeader->Write();
779 }
780
781 ////////////////////////////////////////////////////////////////////////
782
783 void AliUA1JetFinderV1::Init()
784 {
785   // initializes some variables
786
787    // book lego
788   fLego = new
789     TH2F("legoH","eta-phi",
790          fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
791          fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(),
792          fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
793
794 }