Interface to FASTJET. (Rafael Diaz)
[u/mrichter/AliRoot.git] / JETAN / AliFastJetFinder.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 // FastJet finder
18 // interface to FastJet algorithm
19 // Author: Rafael.Diaz.Valdes@cern.ch
20 // kt using NlnN 
21 //---------------------------------------------------------------------
22
23 #include <Riostream.h>
24 #include <TLorentzVector.h>
25 #include <TFile.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TArrayF.h>
29 #include <TRandom.h>
30 #include <TClonesArray.h>
31
32 #include "AliFastJetFinder.h"
33 #include "AliFastJetHeader.h"
34 #include "AliJetReaderHeader.h"
35 #include "AliJetReader.h"
36 #include "AliJet.h"
37 //for FastJet finder
38 #include "FjPseudoJet.hh"
39 #include "FjClusterSequence.hh"
40
41
42 ClassImp(AliFastJetFinder);
43
44
45 ////////////////////////////////////////////////////////////////////////
46
47 AliFastJetFinder::AliFastJetFinder():
48     AliJetFinder(),
49     fHeader(0x0),
50     fLego(0x0),
51     fLegoSignal(0x0)
52 {
53   // Constructor
54 }
55
56 ////////////////////////////////////////////////////////////////////////
57
58 AliFastJetFinder::~AliFastJetFinder()
59
60 {
61   // destructor
62 }
63
64 ////////////////////////////////////////////////////////////////////////
65
66
67 void AliFastJetFinder::FindJets()
68
69 {
70   //create cells and jets array
71   // 1) 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 particles input
77   Float_t* enT  = new Float_t[nIn];
78   Float_t* ptT  = new Float_t[nIn];
79   Float_t* etaT = new Float_t[nIn];
80   Float_t* phiT = new Float_t[nIn];
81   Int_t*   injet = new Int_t[nIn];
82
83
84   //total energy in array
85   Float_t  EtTotal = 0.0;
86   Float_t  meanptCell = 0.0;
87   Float_t  sqptCell = 0.0;
88
89
90   // load input vectors in fLego
91   for (Int_t i = 0; i < nIn; i++){
92     TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
93     enT[i]  = lv->Energy();
94     ptT[i]  = lv->Pt();
95     etaT[i] = lv->Eta();
96     phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
97     if (fReader->GetCutFlag(i) == 1){
98       fLego->Fill(etaT[i], phiT[i], ptT[i]);
99       if(fReader->GetSignalFlag(i) == 1)
100         fLegoSignal->Fill(etaT[i], phiT[i], ptT[i]);
101       EtTotal= EtTotal+ptT[i];
102     }
103   }
104   fJets->SetNinput(nIn);
105
106   // add soft background fixed
107   Int_t nsoft = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
108   Float_t* ptRndm =  new Float_t[nsoft];
109   if(fHeader->AddSoftBackg()){
110     gRandom->RndmArray(nsoft,ptRndm);
111     for(Int_t isoft = 0; isoft < nsoft; isoft++){
112         Float_t ptsoft  = 0.005*ptRndm[isoft];
113         EtTotal= EtTotal+ptsoft;
114     }
115   }
116
117   if(EtTotal == 0){
118      delete enT;
119      delete ptT;
120      delete etaT;
121      delete phiT;
122      return;
123   }
124
125   //cell array
126   Float_t* etCell = new Float_t[90000];   //! Cell Energy // check enough space! *to be done*
127   Float_t* etaCell = new Float_t[90000];  //! Cell eta
128   Float_t* phiCell = new Float_t[90000];  //! Cell phi
129   Int_t*   jetflagCell = new Int_t[90000]; //! Cell flag for jets
130   Float_t* etsigCell = new Float_t[90000];   // signal in this cell
131
132   //jets array
133   Float_t* etaJet = new Float_t[200];
134   Float_t* phiJet = new Float_t[200];
135   Float_t* etJet  = new Float_t[200];
136   Float_t* etsigJet  = new Float_t[200]; //signal energy
137   Float_t* etallJet = new Float_t[200];  //total energy in jet area
138   Int_t*   ncellsJet = new Int_t[200];
139   memset(etaJet,0,sizeof(Float_t)*200);
140   memset(phiJet,0,sizeof(Float_t)*200);
141   memset(etJet,0,sizeof(Float_t)*200);
142   memset(etsigJet,0,sizeof(Float_t)*200);
143   memset(etallJet,0,sizeof(Float_t)*200);
144   memset(ncellsJet,0,sizeof(Int_t)*200);
145
146
147
148   // load cells arrays
149   Int_t nCell = 0;
150   TAxis* xaxis = fLego->GetXaxis();
151   TAxis* yaxis = fLego->GetYaxis();
152   Float_t e = 0.0; Float_t esig = 0.0;
153
154   for (Int_t ib = 1; ib <= fHeader->GetLegoNbinEta(); ib++) {
155       for (Int_t jb = 1; jb <= fHeader->GetLegoNbinPhi(); jb++) {
156                e = fLego->GetBinContent(ib,jb);
157                if (e < 0.0) continue; // don't include this cells
158                Float_t eta  = xaxis->GetBinCenter(ib);
159                Float_t phi  = yaxis->GetBinCenter(jb);
160           if(fHeader->AddSoftBackg())
161              etCell[nCell]  = e + 0.005*ptRndm[nCell];
162           else
163              etCell[nCell]  = e;
164           sqptCell = sqptCell + etCell[nCell]*etCell[nCell]; // xi^2 ////////
165                etaCell[nCell] = eta;
166                phiCell[nCell] = phi;
167           jetflagCell[nCell] = -1; //default
168           esig = fLegoSignal->GetBinContent(ib,jb);
169           if(esig > 0.0)
170              etsigCell[nCell] = esig;
171           else
172              etsigCell[nCell] = 0.0;
173                nCell++;
174       }
175   }
176
177   meanptCell = EtTotal/(Float_t)nCell;
178   sqptCell = sqptCell/(Float_t)nCell;
179
180   Int_t nJets = 0;
181   //call to FastJet Algorithm
182   RunAlgorithm(nJets,etJet,etaJet,phiJet,etsigJet,etallJet,ncellsJet,
183                nCell,etCell,etaCell,phiCell,etsigCell,jetflagCell);
184
185
186   //subtract background
187   SubtractBackg(nCell,jetflagCell,etCell,
188                 nJets,etJet,etallJet,ncellsJet,
189                 meanptCell,sqptCell,EtTotal);
190
191
192   // add jets to list
193   Int_t* index = new Int_t[nJets];
194   Int_t nj = 0;
195   for(Int_t kj=0; kj<nJets; kj++){
196       if ((etaJet[kj] > (fHeader->GetJetEtaMax())) ||
197           (etaJet[kj] < (fHeader->GetJetEtaMin())) ||
198           (etJet[kj] < fHeader->GetMinJetEt())) continue; // acceptance eta range and etmin
199       Float_t px, py,pz,en; // convert to 4-vector
200       px = etJet[kj] * TMath::Cos(phiJet[kj]);
201       py = etJet[kj] * TMath::Sin(phiJet[kj]);
202       pz = etJet[kj] / TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaJet[kj])));
203       en = TMath::Sqrt(px * px + py * py + pz * pz);
204       fJets->AddJet(px, py, pz, en);
205       index[nj] = kj;
206       nj++;
207   }
208
209   //add signal percentage and total signal  in AliJets for analysis tool
210   Float_t* percentage  = new Float_t[nj];
211   Int_t* ncells      = new Int_t[nj];
212   Int_t* mult        = new Int_t[nj];
213
214   for(Int_t i = 0; i< nj; i++){
215      percentage[i] = etsigJet[index[i]]/etJet[index[i]];
216      ncells[i] = ncellsJet[index[i]];
217   }
218
219    //reorder injet flags
220   for(Int_t ipar = 0; ipar < nIn; ipar++){
221      Float_t injetflag =0;
222      Int_t iparCell = fLego->FindBin(etaT[ipar], phiT[ipar]);
223      injet[ipar] = jetflagCell[iparCell];
224      for(Int_t js = 0; js < nj; js++){
225         if(injet[ipar] == index[js]){
226           injet[ipar] = js;  // set the new jet id value
227           mult[js]++; // add multiplicity in jet js
228           injetflag = 1;
229           break;
230         }
231      }
232      if(injetflag == 0) injet[ipar] = -1; // set default value
233   }
234
235
236   fJets->SetNCells(ncells);
237   fJets->SetPtFromSignal(percentage);
238   fJets->SetMultiplicities(mult);
239   fJets->SetInJet(injet);
240   fJets->SetEtaIn(etaT);
241   fJets->SetPhiIn(phiT);
242   fJets->SetPtIn(ptT);
243   fJets->SetEtAvg(meanptCell);
244
245    //delete
246   delete enT;
247   delete ptT;
248   delete etaT;
249   delete phiT;
250   delete injet;
251   //cells
252   delete etCell;
253   delete etaCell;
254   delete phiCell;
255   delete jetflagCell;
256   delete etsigCell;
257   //jets
258   delete etaJet;
259   delete phiJet;
260   delete etJet;
261   delete etsigJet;
262   delete etallJet;
263   delete ncellsJet;
264
265   delete index;
266   delete percentage;
267   delete ncells;
268   delete mult;
269   delete ptRndm;
270
271 }
272
273 ////////////////////////////////////////////////////////////////////////
274 void AliFastJetFinder::RunAlgorithm(Int_t& nJets,Float_t* etJet,Float_t* etaJet,Float_t* phiJet,
275                                     Float_t* etsigJet, Float_t* etallJet, Int_t* ncellsJet,
276                                     Int_t& nCell,Float_t* etCell,Float_t* etaCell,Float_t* phiCell,
277                                     Float_t* etsigCell, Int_t* jetflagCell)
278 {
279    //FastJet objects
280    vector<FjPseudoJet> input_cells; // create a vector
281    for (Int_t i = 0; i < nCell; i++){
282       if(etCell[i] == 0.0) continue; // not include cell empty
283       Double_t px, py,pz,en; // convert to 4-vector
284       px = etCell[i]*TMath::Cos(phiCell[i]);
285       py = etCell[i]*TMath::Sin(phiCell[i]);
286       pz = etCell[i]/TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-etaCell[i])));
287       en = TMath::Sqrt(px * px + py * py + pz * pz);
288       FjPseudoJet input_cell(px,py,pz,en); // create FjPseudoJet object
289       input_cell.set_user_index(i); //label the cell into Fastjet algortihm
290       //push FjPseudoJet of (px,py,pz,en) onto back of the input_cells
291       input_cells.push_back(input_cell);
292    }
293
294    //run the jet clustering with option R=1.0 and strategy= Best
295    Double_t Rparam = fHeader->GetRadius(); // default 1.0;
296    FjClusterSequence clust_seq(input_cells,Rparam);
297
298
299    //vector to get clusters
300    vector<FjPseudoJet> clusters;
301
302    ///////////////////////////////////////////////////////////////////////////
303    //extract the inclusive jets with pt> ptmin sorted by pt
304    //Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
305   // Double_t ptbgRc = EtBackgT*(Rparam*Rparam*TMath::Pi()/areaT);
306   // Double_t ptbgRcfluct = dEtTotal*Rparam*TMath::Sqrt(TMath::Pi()/areaT);
307   // Double_t ptmin = ptbgRc + ptbgRcfluct;
308    clusters = clust_seq.inclusive_jets(0);
309    //////////////////////////////////////////////////////////////////////////
310
311    /////////////////////////////////////////////////////////////////////////
312    //extract the exclusive jets with dcut = 25 GeV**2 and sort them in order
313    //of increasing pt
314    //Float_t areaT = 4*(fHeader->GetLegoEtaMax())*TMath::Pi();
315    //Double_t ptbgRc = EtBackgT*(Rparam*Rparam*TMath::Pi()/areaT);
316    //Double_t ptbgRcfluct = dEtTotal*Rparam*TMath::Sqrt(TMath::Pi()/areaT);
317    //Double_t ptmin = ptbgRc + ptbgRcfluct;
318    //Double_t ktbackg = (fHeader->GetKfactor())*ptmin;
319    //Double_t dcut = ktbackg*ktbackg;
320    //clusters = sorted_by_pt(clust_seq.exclusive_jets(dcut));
321    //clusters = sorted_by_pt(clust_seq.exclusive_jets(5));
322    /////////////////////////////////////////////////////////////////////////
323
324
325    //cout << " ktClusters " << clusters.size() << endl;
326    nJets = (Int_t)clusters.size();
327    ////////////////////////////////////////////////////////////////////////
328    // add all clusters to jet arrays
329    for(UInt_t ij = 0; ij < clusters.size(); ij++){
330        //constituents
331        vector<FjPseudoJet> constituents = clust_seq.constituents(clusters[ij]);
332        //fill jet array info
333        ncellsJet[ij] = (Int_t)constituents.size();
334        phiJet[ij] = clusters[ij].phi();
335        Float_t angle = TMath::ATan(clusters[ij].perp()/clusters[ij].pz());
336        angle = ((angle < 0) ? angle + TMath::Pi() : angle);
337        etaJet[ij] = - TMath::Log(TMath::Tan(angle/2.0));
338        etJet[ij] = clusters[ij].perp();
339        //get constituents cells
340        for(UInt_t jc = 0; jc < constituents.size(); jc++){ // loop for all cells in ij cluster
341            Int_t jcell = constituents[jc].user_index();
342            jetflagCell[jcell] = ij; //flag this cell for jet
343            etsigJet[ij] = etsigJet[ij] + etsigCell[jcell]; // add signal of this cell
344            etallJet[ij] = etallJet[ij] + etCell[jcell];   // add total of this cell
345        }
346    }
347
348 }
349 ////////////////////////////////////////////////////////////////////////
350 void AliFastJetFinder::SubtractBackg(Int_t& nCell, Int_t* jetflagCell, Float_t* etCell,
351                                      Int_t& nJets, Float_t* etJet, Float_t* etallJet, Int_t* ncellsJet,
352                                      Float_t& meanptCell, Float_t& sqptCell, Float_t& etBackg)
353 {
354    // simplest method: subtract background from external region to jets
355
356    //tmp array to flag jets
357    Int_t flagJet[200];
358    //tmp array to flag jet-cell status
359    Int_t tmpjetflagCell[90000];
360
361    Float_t etBackgOld = 0;
362    Float_t prec  = fHeader->GetPrecBg();
363    Float_t bgprec = 1;
364
365    while(bgprec > prec){
366         //clear tmpjetflagCell
367         memset(tmpjetflagCell,-1,sizeof(Int_t)*90000); // init with -1 (all cells are background)
368         //clear flagjet
369         memset(flagJet,0,sizeof(Int_t)*200); // init with 0 (no flag jets)
370         // select clusters > meantmpCell
371         for(Int_t i = 0; i < nJets; i++){
372             Float_t iptcell = etallJet[i]/(Float_t)ncellsJet[i];
373             if(iptcell < meanptCell) continue; // cluster not selected
374             // convert tmp cell background to jet cell
375             for(Int_t ic = 0; ic < nCell; ic++){ //loop for all cells
376                 if(jetflagCell[ic] != i) continue; // other cells
377                 tmpjetflagCell[ic] = i; // convert to a jet cell
378             }
379             //load total energy in cluster
380             etJet[i] = etallJet[i];
381             flagJet[i] = 1; // flag jet
382        }
383        //subtract background
384        for(Int_t j = 0; j < nCell; j++){ // loop for all cells
385            Int_t idxjet = tmpjetflagCell[j];
386            if(idxjet == -1) continue; // background cell
387            if(idxjet == -2) continue; // background protected cell
388            etJet[idxjet] = etJet[idxjet] - meanptCell;
389        }
390        // evaluate background fluctuations (rms value)
391        Float_t rmsptCell = TMath::Sqrt(sqptCell - meanptCell*meanptCell);
392        //fake jets
393        for(Int_t k = 0; k < nJets; k++){
394            if(flagJet[k] != 1) continue; // only flaged jets
395            //if(etJet[k] > fHeader->GetMinJetEt()) continue;  // jet ok!!
396            if(etJet[k] > rmsptCell*ncellsJet[k]) continue;  // jet ok!!
397            //clear tmpjetflag in cells
398            for(Int_t kc = 0; kc < nCell; kc++){ //loop for all cells
399                if(tmpjetflagCell[kc] != k) continue; // other cells
400                tmpjetflagCell[kc] = -1; // convert to background tmp cell
401            }
402            // clear all previous jet flags
403            etJet[k] = 0;
404            flagJet[k] = 0;
405        }
406        // recalculate background
407        sqptCell = 0;
408        etBackgOld = etBackg;
409        etBackg = 0;
410        Int_t nCellBackg = 0;
411        for(Int_t l = 0; l < nCell; l++){ // loop for all cells
412           if(tmpjetflagCell[l] != -1) continue; // cell included in some jet or protected
413           nCellBackg++;
414           etBackg = etBackg + etCell[l]; //add cell to background
415           //calc sqptCell
416           sqptCell = sqptCell + etCell[l]*etCell[l];
417        }
418        if(nCellBackg){
419           meanptCell = etBackg/(Float_t)nCellBackg; // new pt cell mean value
420           sqptCell = sqptCell/(Float_t)nCellBackg;
421        }else{
422           meanptCell = 0;
423           sqptCell = 0;
424        }
425        // evaluate presicion values
426        if(etBackg)
427          bgprec = (etBackgOld - etBackg)/etBackg;
428        else
429          bgprec = 0;
430    }
431
432    // set etJet 0 for all clusters not flaged in order to
433    for(Int_t m = 0; m < nJets; m++){
434        if(flagJet[m] == 1) continue; // flaged jets
435        etJet[m] = 0; //others clusters
436    }
437
438
439 }
440
441 ////////////////////////////////////////////////////////////////////////
442 void AliFastJetFinder::SubtractBackgArea(Int_t& nCell, Int_t* jetflagCell,
443                                      Float_t* etCell,Int_t&nJets, Float_t* etJet, Float_t* etallJet)
444 {
445    // area method: subtract background from external region to jets
446    // using fixed area pi*Rc2
447
448    // n cells contained in a cone Rc
449    Double_t Rc = fHeader->GetRadius();
450    Float_t nCellRc = Rc*Rc*TMath::Pi()/(0.015*0.015); // change in future !!!!
451    //tmp array to flag fake jets
452    Int_t fakeflagJet[100];
453    memset(fakeflagJet,0,sizeof(Int_t)*100); // init with 0 (no fake jets)
454    Int_t njfake = nJets;
455    while(njfake > 0){
456        //calc background per cell
457        Int_t nCellBackg = 0;
458        Float_t EtBackg = 0.0;
459        for(Int_t i = 0; i < nCell; i++){ // loop for all cells
460            if(jetflagCell[i] != -1) continue; // cell included in some jet
461            nCellBackg++;
462            EtBackg = EtBackg + etCell[i]; //add cell to background
463        }
464        //subtract background energy per jet
465        for(Int_t l = 0; l < nJets; l++){
466            if(fakeflagJet[l] == 1) continue; // fake jet
467            etJet[l] = etallJet[l] - nCellRc*EtBackg/(Float_t)nCellBackg;
468        }
469        //fake jets analysis
470        njfake = 0;
471        for(Int_t k = 0; k < nJets; k++){
472            if(fakeflagJet[k] == 1) continue; // fake jet
473            if(etJet[k] < fHeader->GetMinJetEt()){  // a new fake jet
474               //clear jet flag in cells
475               for(Int_t kc = 0; kc < nCell; kc++){ //loop for all cells
476                   Int_t kidx = jetflagCell[kc];
477                   if(kidx != k) continue; // other cells
478                   jetflagCell[kc] = -1; // convert to background cell
479               }
480               fakeflagJet[k] = 1; // mark as a fake jet
481               njfake++; //count fakes in this loop
482            }
483        }
484    }
485 }
486
487 ////////////////////////////////////////////////////////////////////////
488
489
490
491 void AliFastJetFinder::Reset()
492 {
493   fLego->Reset();
494   fLegoSignal->Reset();
495   fJets->ClearJets();
496
497 }
498 ////////////////////////////////////////////////////////////////////////
499
500 void AliFastJetFinder::WriteJHeaderToFile()
501 {
502   fOut->cd();
503   fHeader->Write();
504 }
505
506 ////////////////////////////////////////////////////////////////////////
507
508 void AliFastJetFinder::Init()
509 {
510   // initializes some variables
511    // book lego
512   fLego = new
513     TH2F("legoH","eta-phi",
514          fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
515          fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(),
516          fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
517   fLegoSignal = new
518     TH2F("legoSignalH","eta-phi signal",
519          fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
520          fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(),
521          fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
522
523 }