723417d36c8edb20f305f7f5bf442b4acec0a98e
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinder.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 jet studies
20 // manages the search for jets using charged particle momentum and 
21 // neutral cell energy information
22 // Authors: Rafael.Diaz.Valdes@cern.ch 
23 //          magali.estienne@subatech.in2p3.fr
24 //          alexandre.shabetai@cern.ch
25 // ** 2011
26 // Modified accordingly to reader/finder splitting and new handling of neutral information
27 // Versions V1 and V2 merged
28 //---------------------------------------------------------------------
29
30 #include <TH2F.h>
31 #include <TMath.h>
32
33 #include "AliUA1JetFinder.h"
34 #include "AliUA1JetHeaderV1.h"
35 #include "AliJetCalTrk.h"
36 #include "AliJetBkg.h"
37 #include "AliAODJetEventBackground.h"
38 #include "AliAODJet.h"
39
40 ClassImp(AliUA1JetFinder)
41
42 ////////////////////////////////////////////////////////////////////////
43
44 AliUA1JetFinder::AliUA1JetFinder():
45   AliJetFinder(),
46   fLego(0),  
47   fJetBkg(new AliJetBkg())
48 {
49   // Default constructor
50 }
51
52 //-----------------------------------------------------------------------
53 AliUA1JetFinder::~AliUA1JetFinder()
54 {
55   // Destructor
56   delete fLego;
57   delete fJetBkg;
58
59 }
60
61 //-----------------------------------------------------------------------
62 void AliUA1JetFinder::FindJets()
63 {
64   //  Used to find jets using charged particle momentum information 
65   //  & neutral energy from calo cells
66   //
67   //  1) Fill cell map array
68   //  2) calculate total energy and fluctuation level
69   //  3) Run algorithm
70   //     3.1) look centroides in cell map
71   //     3.2) calculate total energy in cones
72   //     3.3) flag as a possible jet
73   //     3.4) reorder cones by energy
74   //  4) subtract backg in accepted jets
75   //  5) fill AliJet list
76   
77   // transform input to pt,eta,phi plus lego
78   
79   AliUA1JetHeader* header  = (AliUA1JetHeader*) fHeader;
80   Int_t              nIn = fCalTrkEvent->GetNCalTrkTracks();
81   fDebug   = fHeader->GetDebug();
82   
83   if (nIn <= 0) return;
84   fJetBkg->SetHeader(fHeader);
85   fJetBkg->SetCalTrkEvent(GetCalTrkEvent());
86   fJetBkg->SetDebug(fDebug);
87   // local arrays for input
88   // ToDo: check memory fragmentation, maybe better to
89   // define them globally and resize as needed
90   // Fragmentation should be worse for low mult...
91   Float_t* ptT      = new Float_t[nIn];
92   Float_t* etaT     = new Float_t[nIn];
93   Float_t* phiT     = new Float_t[nIn];
94   Int_t*   injet    = new Int_t[nIn];
95   Int_t*   injetOk  = new Int_t[nIn];
96
97   memset(ptT,0,sizeof(Float_t)*nIn);
98   memset(etaT,0,sizeof(Float_t)*nIn);
99   memset(phiT,0,sizeof(Float_t)*nIn);
100   memset(injet,0,sizeof(Int_t)*nIn);
101   memset(injetOk,-1,sizeof(Int_t)*nIn);
102
103   // load input vectors and calculate total energy in array
104
105   // total energy in array
106   Float_t etbgTotal = 0.;
107   Float_t npart = 0.;
108   Float_t etbg2 = 0.;
109
110   for (Int_t i = 0; i < fCalTrkEvent->GetNCalTrkTracks(); i++){
111     ptT[i]  = fCalTrkEvent->GetCalTrkTrack(i)->GetPt();
112     etaT[i] = fCalTrkEvent->GetCalTrkTrack(i)->GetEta();
113     phiT[i] = ((fCalTrkEvent->GetCalTrkTrack(i)->GetPhi() < 0) ? (fCalTrkEvent->GetCalTrkTrack(i)->GetPhi()) + 2 * TMath::Pi() : fCalTrkEvent->GetCalTrkTrack(i)->GetPhi());
114     //fCalTrkEvent->GetCalTrkTrack(i)->Print(Form("%d",i));
115     if (fCalTrkEvent->GetCalTrkTrack(i)->GetCutFlag() != 1) continue;
116     fLego->Fill(etaT[i], phiT[i], ptT[i]);
117     npart += 1;
118     etbgTotal+= ptT[i];
119     etbg2 += ptT[i]*ptT[i];
120   }
121   
122   // calculate total energy and fluctuation in map
123   Double_t meanpt = 0.;
124   Double_t ptRMS = 0.;
125   if(npart>0){
126     meanpt = etbgTotal/npart;
127     etbg2 = etbg2/npart;
128     if(etbg2>(meanpt*meanpt)){// prenent NAN, should only happen due to numerical instabilities
129       ptRMS = TMath::Sqrt(etbg2-meanpt*meanpt);
130     }
131   }
132   Double_t dEtTotal = (TMath::Sqrt(npart))*TMath::Sqrt(meanpt * meanpt + ptRMS*ptRMS);
133   
134   // arrays to hold jets
135   Float_t etaJet[kMaxJets];
136   Float_t phiJet[kMaxJets];
137   Float_t etJet[kMaxJets];
138   Float_t etsigJet[kMaxJets]; //signal et in jet
139   Float_t etallJet[kMaxJets];  // total et in jet (tmp variable)
140   Int_t   ncellsJet[kMaxJets];
141   Int_t   multJetT[kMaxJets];
142   Int_t   multJetC[kMaxJets];
143   Int_t   multJet[kMaxJets];
144   Float_t *areaJet = new Float_t[kMaxJets];
145   // Used for jet reordering at the end of the jet finding procedure
146   Float_t etaJetOk[kMaxJets];
147   Float_t phiJetOk[kMaxJets];
148   Float_t etJetOk[kMaxJets];
149   Float_t etsigJetOk[kMaxJets]; //signal et in jet
150   Float_t etallJetOk[kMaxJets];  // total et in jet (tmp variable)
151   Int_t   ncellsJetOk[kMaxJets];
152   Int_t   multJetOk[kMaxJets];
153   Float_t *areaJetOk = new Float_t[kMaxJets];
154   Int_t   nJets; // to hold number of jets found by algorithm
155   Int_t   nj;    // number of jets accepted
156   Float_t prec  = header->GetPrecBg();
157   Float_t bgprec = 1;
158
159   while(bgprec > prec){
160     //reset jet arrays in memory
161     memset(etaJet,0,sizeof(Float_t)*kMaxJets);
162     memset(phiJet,0,sizeof(Float_t)*kMaxJets);
163     memset(etJet,0,sizeof(Float_t)*kMaxJets);
164     memset(etallJet,0,sizeof(Float_t)*kMaxJets);
165     memset(etsigJet,0,sizeof(Float_t)*kMaxJets);
166     memset(ncellsJet,0,sizeof(Int_t)*kMaxJets);
167     memset(multJetT,0,sizeof(Int_t)*kMaxJets);
168     memset(multJetC,0,sizeof(Int_t)*kMaxJets);
169     memset(multJet,0,sizeof(Int_t)*kMaxJets);
170     memset(areaJet,0,sizeof(Float_t)*kMaxJets);
171     memset(etaJetOk,0,sizeof(Float_t)*kMaxJets);
172     memset(phiJetOk,0,sizeof(Float_t)*kMaxJets);
173     memset(etJetOk,0,sizeof(Float_t)*kMaxJets);
174     memset(etallJetOk,0,sizeof(Float_t)*kMaxJets);
175     memset(etsigJetOk,0,sizeof(Float_t)*kMaxJets);
176     memset(ncellsJetOk,0,sizeof(Int_t)*kMaxJets);
177     memset(multJetOk,0,sizeof(Int_t)*kMaxJets);
178     memset(areaJetOk,0,sizeof(Float_t)*kMaxJets);
179     nJets = 0;
180     nj = 0;
181     // reset particles-jet array in memory
182     memset(injet,-1,sizeof(Int_t)*nIn);
183     //run cone algorithm finder
184     RunAlgoritm(etbgTotal,dEtTotal,nJets,etJet,etaJet,phiJet,etallJet,ncellsJet);
185     //run background subtraction
186     if(nJets > header->GetNAcceptJets()) // limited number of accepted jets per event
187       nj = header->GetNAcceptJets();
188     else
189       nj = nJets;
190
191     //subtract background
192     Float_t etbgTotalN = 0.0; //new background
193     Float_t sigmaN = 0.0; //new background
194     if(header->GetBackgMode() == 1) {// standard
195       fJetBkg->SubtractBackg(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);}
196     if(header->GetBackgMode() == 2) //cone
197       fJetBkg->SubtractBackgCone(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
198     if(header->GetBackgMode() == 3) //ratio
199       fJetBkg->SubtractBackgRatio(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
200     if(header->GetBackgMode() == 4) //statistic
201       fJetBkg->SubtractBackgStat(nIn,nj,etbgTotalN,sigmaN,ptT,etaT,phiT,etJet,etaJet,phiJet,etsigJet,multJetT,multJetC,multJet,injet,areaJet);
202
203     //calc precision
204     if(etbgTotalN != 0.0)
205       bgprec = (etbgTotal - etbgTotalN)/etbgTotalN;
206     else
207       bgprec = 0;
208     etbgTotal = etbgTotalN; // update with new background estimation
209
210   } //end while
211
212   // add tracks to the jet if it wasn't yet done
213   if (header->GetBackgMode() == 0){
214     Float_t rc= header->GetRadius();
215     for(Int_t jpart = 0; jpart < nIn; jpart++){ // loop for all particles in array
216       for(Int_t ijet=0; ijet<nJets; ijet++){
217         Float_t deta = etaT[jpart] - etaJet[ijet];
218         Float_t dphi = phiT[jpart] - phiJet[ijet];
219         if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
220         if (dphi >  TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
221         Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
222         if(dr <= rc){ // particles inside this cone
223           injet[jpart] = ijet;
224           break;
225         }
226       }// end jets loop
227     } //end particle loop
228   }
229
230   // add jets to list
231   if (fDebug>1) printf("Found %d jets \n", nj);
232   
233   // Reorder jets by et in cone
234   // Sort jets by energy
235   Int_t idx[kMaxJets];
236   TMath::Sort(nJets, etJet, idx);
237   for(Int_t p = 0; p < nJets; p++)
238     {
239       etaJetOk[p]    = etaJet[idx[p]];
240       phiJetOk[p]    = phiJet[idx[p]];
241       etJetOk[p]     = etJet[idx[p]];
242       etallJetOk[p]  = etJet[idx[p]];
243       etsigJetOk[p]  = etsigJet[idx[p]];
244       ncellsJetOk[p] = ncellsJet[idx[p]];
245       multJetOk[p]   = multJet[idx[p]];
246       areaJetOk[p]   = areaJet[idx[p]];
247     }
248
249   //////////////////////////
250
251   Int_t nTracks = fCalTrkEvent->GetNCalTrkTracks();
252   
253   for(Int_t kj=0; kj<nj; kj++)
254     {
255       if ((etaJetOk[kj] > (header->GetJetEtaMax())) ||
256           (etaJetOk[kj] < (header->GetJetEtaMin())) ||
257           (phiJetOk[kj] > (header->GetJetPhiMax())) ||
258           (phiJetOk[kj] < (header->GetJetPhiMin())) ||
259           (etJetOk[kj] < header->GetMinJetEt())) continue; // acceptance eta range and etmin
260       Float_t px=-999, py=-999 ,pz=-999 ,en=-999; // convert to 4-vector
261       px = etJetOk[kj] * TMath::Cos(phiJetOk[kj]);
262       py = etJetOk[kj] * TMath::Sin(phiJetOk[kj]);
263       pz = etJetOk[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJetOk[kj])));
264       en = TMath::Sqrt(px * px + py * py + pz * pz);
265       AliAODJet jet(px, py, pz, en);
266
267       // Calc jet area if it wasn't yet done
268       if (header->GetBackgMode() == 0){
269         // calculate the area of the jet
270         Float_t rc= header->GetRadius();
271         areaJetOk[kj] = fJetBkg->CalcJetAreaEtaCut(rc,etaJetOk[kj]);
272       }
273
274       jet.SetEffArea(areaJetOk[kj],0.,0.,0.);
275       jet.SetBgEnergy(etbgTotal,0.);
276       if (fDebug>1) jet.Print(Form("%d",kj));
277       
278       for(Int_t jpart = 0; jpart < nTracks; jpart++) { // loop for all particles in array
279         // Track to jet reordering
280         if(injet[jpart] == idx[kj]){
281           injetOk[jpart] = kj;
282         }
283         // Check if the particle belongs to the jet and add the ref
284         if(injetOk[jpart] == kj && fCalTrkEvent->GetCalTrkTrack(jpart)->GetCutFlag() == 1) {
285           jet.AddTrack(fCalTrkEvent->GetCalTrkTrack(jpart)->GetTrackObject());
286         }
287       }
288
289       AddJet(jet);
290       
291     }
292
293   //delete
294   delete[] ptT;
295   delete[] etaT;
296   delete[] phiT;
297   delete[] injet;
298   delete[] injetOk;
299   delete[] areaJet;
300   delete[] areaJetOk;
301
302 }
303
304 //-----------------------------------------------------------------------
305 void AliUA1JetFinder::RunAlgoritm(Float_t etbgTotal, Double_t dEtTotal, Int_t& nJets,
306                                   Float_t* const etJet,Float_t* const etaJet, Float_t* const phiJet,
307                                   Float_t* const etallJet, Int_t* const ncellsJet)
308 {
309   // Dump lego
310   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
311   const Int_t nBinsMax = 120000; // we use a fixed array not to fragment memory 
312  
313   const Int_t nBinEta = header->GetLegoNbinEta();
314   const Int_t nBinPhi = header->GetLegoNbinPhi();
315   if((nBinPhi*nBinEta)>nBinsMax){
316     AliError("Too many bins of the ETA-PHI histogram");
317   }
318
319   Float_t etCell[nBinsMax] = {0.};   //! Cell Energy
320   Float_t etaCell[nBinsMax] = {0.};  //! Cell eta
321   Float_t phiCell[nBinsMax] = {0.};  //! Cell phi
322   Short_t flagCell[nBinsMax] = {0}; //! Cell flag
323   
324   Int_t nCell = 0;
325   TAxis* xaxis = fLego->GetXaxis();
326   TAxis* yaxis = fLego->GetYaxis();
327   Float_t e = 0.0;
328   for (Int_t i = 1; i <= nBinEta; i++) {
329     for (Int_t j = 1; j <= nBinPhi; j++) {
330       e = fLego->GetBinContent(i,j);
331       if (e < 0.0) continue; // don't include this cells
332       Float_t eta  = xaxis->GetBinCenter(i);
333       Float_t phi  = yaxis->GetBinCenter(j);
334       etCell[nCell]  = e;
335       etaCell[nCell] = eta;
336       phiCell[nCell] = phi;
337       flagCell[nCell] = 0; //default
338       nCell++;
339     }
340   }
341   // Parameters from header
342   Float_t minmove = header->GetMinMove();
343   Float_t maxmove = header->GetMaxMove();
344   Float_t rc      = header->GetRadius();
345   Float_t etseed  = header->GetEtSeed();
346   // Tmp array of jets form algoritm
347   Float_t etaAlgoJet[kMaxJets] = {0.0};
348   Float_t phiAlgoJet[kMaxJets] = {0.0};
349   Float_t etAlgoJet[kMaxJets] = {0.0};
350   Int_t   ncellsAlgoJet[kMaxJets] = {0};
351
352   // Run algorithm//
353   
354   // Sort cells by et
355   Int_t index[nBinsMax];
356   TMath::Sort(nCell, etCell, index);
357   // variable used in centroide loop
358   Float_t eta   = 0.0;
359   Float_t phi   = 0.0;
360   Float_t eta0  = 0.0;
361   Float_t phi0  = 0.0;
362   Float_t etab  = 0.0;
363   Float_t phib  = 0.0;
364   Float_t etas  = 0.0;
365   Float_t phis  = 0.0;
366   Float_t ets   = 0.0;
367   Float_t deta  = 0.0;
368   Float_t dphi  = 0.0;
369   Float_t dr    = 0.0;
370   Float_t etsb  = 0.0;
371   Float_t etasb = 0.0;
372   Float_t phisb = 0.0;
373   Float_t dphib = 0.0;
374
375   for(Int_t icell = 0; icell < nCell; icell++)
376     {
377       Int_t jcell = index[icell];
378       if(etCell[jcell] <= etseed) continue; // if cell energy is low et seed
379       if(flagCell[jcell] != 0) continue; // if cell was used before
380       
381       eta  = etaCell[jcell];
382       phi  = phiCell[jcell];
383       eta0 = eta;
384       phi0 = phi;
385       etab = eta;
386       phib = phi;
387       ets  = etCell[jcell];
388       etas = 0.0;
389       phis = 0.0;
390       etsb = ets;
391       etasb = 0.0;
392       phisb = 0.0;
393       for(Int_t kcell =0; kcell < nCell; kcell++)
394         {
395           Int_t lcell = index[kcell];
396           if(lcell == jcell) continue; // cell itself
397           if(flagCell[lcell] != 0) continue; // cell used before
398           if(etCell[lcell] > etCell[jcell]) continue; // can this happen
399           //calculate dr
400           deta = etaCell[lcell] - eta;
401           dphi = TMath::Abs(phiCell[lcell] - phi);
402           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
403           dr = TMath::Sqrt(deta * deta + dphi * dphi);
404           if(dr <= rc)
405             {
406               // calculate offset from initiate cell
407               deta = etaCell[lcell] - eta0;
408               dphi = phiCell[lcell] - phi0;
409               if (dphi < -TMath::Pi()) dphi= dphi + 2.0 * TMath::Pi();
410               if (dphi > TMath::Pi()) dphi = dphi - 2.0 * TMath::Pi();
411               etas = etas + etCell[lcell]*deta;
412               phis = phis + etCell[lcell]*dphi;
413               ets = ets + etCell[lcell];
414               //new weighted eta and phi including this cell
415               eta = eta0 + etas/ets;
416               phi = phi0 + phis/ets;
417               // if cone does not move much, just go to next step
418               dphib = TMath::Abs(phi - phib);
419               if (dphib > TMath::Pi()) dphib = 2. * TMath::Pi() - dphib;
420               dr = TMath::Sqrt((eta-etab)*(eta-etab) + dphib * dphib);
421               if(dr <= minmove) break;
422               // cone should not move more than max_mov
423               dr = TMath::Sqrt((etas/ets)*(etas/ets) + (phis/ets)*(phis/ets));
424               if(dr > maxmove){
425                 eta = etab;
426                 phi = phib;
427                 ets = etsb;
428                 etas = etasb;
429                 phis = phisb;
430               } else { // store this loop information
431                 etab=eta;
432                 phib=phi;
433                 etsb = ets;
434                 etasb = etas;
435                 phisb = phis;
436               }
437             } // inside cone
438         }//end of cells loop looking centroide
439       
440       // Avoid cones overloap (to be implemented in the future)
441       
442       // Flag cells in Rc, estimate total energy in cone
443       Float_t etCone   = 0.0;
444       Int_t   nCellIn  = 0;
445       rc = header->GetRadius();
446
447       for(Int_t ncell =0; ncell < nCell; ncell++)
448         {
449           if(flagCell[ncell] != 0) continue; // cell used before
450           //calculate dr
451           deta = etaCell[ncell] - eta;
452           dphi = phiCell[ncell] - phi;
453           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
454           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
455           dr = TMath::Sqrt(deta * deta + dphi * dphi);
456           if(dr <= rc){  // cell in cone
457             flagCell[ncell] = -1;
458             etCone+=etCell[ncell];
459             nCellIn++;
460           }
461         }
462       
463       // Select jets with et > background
464       // estimate max fluctuation of background in cone
465       Double_t ncellin = (Double_t)nCellIn;
466       Double_t ntcell  = (Double_t)nCell;
467       Double_t etbmax = (etbgTotal + dEtTotal )*(ncellin/ntcell);
468       // min cone et
469       Double_t etcmin = etCone ;  // could be used etCone - etmin !!
470       //decisions !! etbmax < etcmin
471       
472       for(Int_t mcell =0; mcell < nCell; mcell++){
473         if(flagCell[mcell] == -1){
474           if(etbmax < etcmin)
475             flagCell[mcell] = 1; //flag cell as used
476           else
477             flagCell[mcell] = 0; // leave it free
478         }
479       }
480       //store tmp jet info !!!
481       if(etbmax < etcmin) {
482         if(nJets<kMaxJets){
483           etaAlgoJet[nJets]    = eta;
484           phiAlgoJet[nJets]    = phi;
485           etAlgoJet[nJets]     = etCone;
486           ncellsAlgoJet[nJets] = nCellIn;
487           nJets++;
488         }
489         else{
490           AliError(Form("Too many jets (> %d) found by UA1JetFinder, adapt your cuts",kMaxJets));
491           break;
492         }
493       }      
494     } // end of cells loop
495
496   //reorder jets by et in cone
497   //sort jets by energy
498   for(Int_t p = 0; p < nJets; p++)
499     {
500       etaJet[p]    = etaAlgoJet[p];
501       phiJet[p]    = phiAlgoJet[p];
502       etJet[p]     = etAlgoJet[p];
503       etallJet[p]  = etAlgoJet[p];
504       ncellsJet[p] = ncellsAlgoJet[p];
505     }
506
507 }
508
509 //-----------------------------------------------------------------------
510 void AliUA1JetFinder::Reset()
511 {
512   fLego->Reset();
513   AliJetFinder::Reset();
514
515 }
516
517 //-----------------------------------------------------------------------
518 void AliUA1JetFinder::WriteJHeaderToFile() const
519 {
520   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
521   header->Write();
522
523 }
524
525 //-----------------------------------------------------------------------
526 void AliUA1JetFinder::Init()
527 {
528
529   // initializes some variables
530   AliUA1JetHeader* header = (AliUA1JetHeader*) fHeader;
531   // book lego
532   fLego = new TH2F("legoH","eta-phi",
533                    header->GetLegoNbinEta(), header->GetLegoEtaMin(),
534                    header->GetLegoEtaMax(),  header->GetLegoNbinPhi(),
535                    header->GetLegoPhiMin(),  header->GetLegoPhiMax());
536   
537 }
538