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 }