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