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 //---------------------------------------------------------------------
  17 // UA1 Jet finder
  18 // manages the search for jets
  19 // Author: jgcn@mda.cinvestav.mx
  20 // (code adapted from EMCAL directory)
  21 //---------------------------------------------------------------------
  22
  23 #include <Riostream.h>
  24 #include <TLorentzVector.h>
  25 #include <TFile.h>
  26 #include <TH2F.h>
  27 #include "AliUA1JetFinder.h"
  28 #include "AliUA1JetHeader.h"
  29 #include "UA1Common.h"
  30 #include "AliJetReaderHeader.h"
  31 #include "AliJetReader.h"
  32 #include "AliJet.h"
  33
  34
  35 ClassImp(AliUA1JetFinder)
  36
  37 ////////////////////////////////////////////////////////////////////////
  38
  39 AliUA1JetFinder::AliUA1JetFinder():
  40   fHeader(0x0),
  41   fLego(0x0)
  42 {
  43   // Default constructor
  44 }
  45
  46 ////////////////////////////////////////////////////////////////////////
  47
  48 AliUA1JetFinder::~AliUA1JetFinder()
  49
  50 {
  51   // destructor
  52 }
  53
  54 ////////////////////////////////////////////////////////////////////////
  55
  56 #ifndef WIN32
  57 # define ua1_jet_finder ua1_jet_finder_
  58 # define hf1 hf1_
  59 # define type_of_call
  60
  61 #else
  62 # define ua1_jet_finder UA1_JET_FINDER
  63 # define hf1 HF1
  64 # define type_of_call _stdcall
  65 #endif
  66
  67 extern "C" void type_of_call
  68 ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
  69                Float_t  etc[60000],  Float_t etac[60000],
  70                Float_t  phic[60000],
  71                Float_t& min_move, Float_t& max_move, Int_t& mode,
  72                Float_t& prec_bg,  Int_t& ierror);
  73
  74 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
  75
  76 ////////////////////////////////////////////////////////////////////////
  77
  78 void AliUA1JetFinder::FindJetsTPC()
  79
  80 {
  81   // initialize event, look for jets, download jet info
  82
  83   // initialization in 2 steps
  84   // 1) transform input to pt,eta,phi plus lego
  85   // 2) dump lego
  86
  87   // 1) transform input to pt,eta,phi plus lego
  88   TClonesArray *lvArray = fReader->GetMomentumArray();
  89   Int_t nIn =  lvArray->GetEntries();
  90   if (nIn == 0) return;
  91
  92   // local arrays for input
  93   Float_t* enT  = new Float_t[nIn];
  94   Float_t* ptT  = new Float_t[nIn];
  95   Float_t* etaT = new Float_t[nIn];
  96   Float_t* phiT = new Float_t[nIn];
  97
  98   // load input vectors
  99   for (Int_t i = 0; i < nIn; i++){
 100     TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
 101     enT[i]  = lv->Energy();
 102     ptT[i]  = lv->Pt();
 103     etaT[i] = lv->Eta();
 104     phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
 105     if (fReader->GetCutFlag(i) == 1)
 106       fLego->Fill(etaT[i], phiT[i], ptT[i]);
 107   }
 108   fJets->SetNinput(nIn);
 109
 110   // 2) dump lego
 111   // check enough space! *to be done*
 112   Float_t etCell[60000];   //! Cell Energy
 113   Float_t etaCell[60000];  //! Cell eta
 114   Float_t phiCell[60000];  //! Cell phi
 115
 116   Int_t nCell = 0;
 117   TAxis* xaxis = fLego->GetXaxis();
 118   TAxis* yaxis = fLego->GetYaxis();
 119   Float_t e = 0.0;
 120   for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
 121       for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
 122           e = fLego->GetBinContent(i,j);
 123           if (e > 0.0) e -= fHeader->GetMinCellEt();
 124           if (e < 0.0) e = 0.;
 125           Float_t eta  = xaxis->GetBinCenter(i);
 126           Float_t phi  = yaxis->GetBinCenter(j);
 127           etCell[nCell]  = e;
 128           etaCell[nCell] = eta;
 129           phiCell[nCell] = phi;
 130           nCell++;
 131       }
 132   }
 133
 134   // run the algo. Parameters from header
 135   Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
 136   Float_t minmove = fHeader->GetMinMove();
 137   Float_t maxmove = fHeader->GetMaxMove();
 138   Int_t mode      = fHeader->GetMode();
 139   Float_t precbg  = fHeader->GetPrecBg();
 140   Int_t ierror;
 141   ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell,
 142                  minmove, maxmove, mode, precbg, ierror);
 143
 144   // sort jets
 145   Int_t * idx  = new Int_t[UA1JETS.njet];
 146   TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
 147
 148   // download jet info.
 149   for(Int_t i = 0; i < UA1JETS.njet; i++) {
 150     // reject events outside acceptable eta range
 151     if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
 152         || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
 153           continue;
 154       Float_t px, py,pz,en; // convert to 4-vector
 155       px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
 156       py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
 157       pz = UA1JETS.etj[idx[i]] /
 158           TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
 159       en = TMath::Sqrt(px * px + py * py + pz * pz);
 160       fJets->AddJet(px, py, pz, en);
 161   }
 162
 163   // find multiplicities and relationship jet-particle
 164   // find also percentage of pt from pythia
 165   Int_t* injet = new Int_t[nIn];
 166   Int_t* sflag = new Int_t[nIn];
 167   for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
 168   Int_t* mult  = new Int_t[fJets->GetNJets()];
 169   Int_t* ncell  = new Int_t[fJets->GetNJets()];
 170   Float_t* percentage  = new Float_t[fJets->GetNJets()];
 171
 172   for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
 173       Float_t pt_sig = 0.0;
 174       mult[i] = 0;
 175       ncell[i] = UA1JETS.ncellj[i];
 176       for (Int_t j = 0; j < nIn; j++) {
 177           Float_t deta = etaT[j] - fJets->GetEta(i);
 178           Float_t dphi = phiT[j] - fJets->GetPhi(i);
 179           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 180           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 181           Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 182           if (dr < fHeader->GetRadius() && injet[j] == 0) {
 183               injet[j] = -(i+1);
 184               if ((fReader->GetCutFlag(j)) == 1 &&
 185                   (etaT[j] < fHeader->GetLegoEtaMax()) &&
 186                   (etaT[j] > fHeader->GetLegoEtaMin())) {
 187                   injet[j] = i+1;
 188                   mult[i]++;
 189                   if (fReader->GetSignalFlag(j) == 1) {
 190                     pt_sig+=ptT[j];
 191                     sflag[j]=1;
 192                   }
 193               }
 194           }
 195       }
 196       percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
 197         ((Double_t) fJets->GetPt(i));
 198   }
 199   fJets->SetNCells(ncell);
 200   fJets->SetPtFromSignal(percentage);
 201   fJets->SetMultiplicities(mult);
 202   fJets->SetInJet(injet);
 203   fJets->SetEtaIn(etaT);
 204   fJets->SetPhiIn(phiT);
 205   fJets->SetPtIn(ptT);
 206   fJets->SetEtAvg(UA1JETS.etavg);
 207   delete ncell;
 208   delete enT;
 209   delete ptT;
 210   delete etaT;
 211   delete phiT;
 212   delete injet;
 213   delete idx;
 214   delete mult;
 215   delete percentage;
 216 }
 217
 218 void AliUA1JetFinder::FindJets()
 219 { // Find jets with TPC or EMCAL or TPC + EMCAL information
 220   // initialize event, look for jets, download jet info
 221
 222   // 1) transform input to pt,eta,phi plus lego
 223     AliJetUnitArray *fUnit = fReader->GetUnitArray();
 224     Int_t nIn = fUnit->GetUnitEntries();
 225     Int_t fOption = fReader->GetReaderHeader()->GetDetector();
 226     Int_t fDebug = fReader->GetReaderHeader()->GetDebug();
 227
 228     if(fDebug>1){
 229         printf("Inside FindJets function ! \n");
 230         printf("fOption : %d \n",fOption);
 231     }
 232
 233     if(fDebug>10){
 234         cout << "fUnit : " << fUnit << endl;
 235         printf("nIn = fUnit->GetUnitEntries() : %d \n",nIn);
 236     }
 237
 238     if(fDebug>10){
 239         printf("     -----------------------------------------------------------------\n");
 240         printf("     --- All inputs in fUnitArray ---\n");
 241         for(Int_t i=0; i<nIn ; i++){
 242             if(fUnit[i].GetUnitEnergy()!=0){
 243                 printf("     -----------------------------------------------------------------\n");
 244                 cout << "     |  ID : " << fUnit[i].GetUnitID() << "  |  Eta : " << fUnit[i].GetUnitEta() << "  |  Phi : " << fUnit[i].GetUnitPhi() << "  |  Energy : " << fUnit[i].GetUnitEnergy() << " |" <<endl;
 245             }
 246         }
 247     }
 248
 249     if (nIn == 0) return;
 250
 251     Int_t nCandidateCut = 0;
 252     Int_t nCandidateCut2 = 0;
 253     Int_t nCandidate = 0;
 254     for (Int_t i = 0; i < nIn; i++){
 255         if(fUnit[i].GetUnitEnergy()>0. && fUnit[i].GetUnitCutFlag()==1) {
 256             // Number of good tracks/digits which have passed the pt cut
 257             nCandidateCut += 1;
 258         }
 259     if(fUnit[i].GetUnitEnergy()>0.) {
 260         // Number of good tracks/digits without pt cut
 261         nCandidate += 1;
 262     }
 263     }
 264
 265     if(fDebug>=3){
 266     cout << "nCandidate : " << nCandidate << endl;
 267     cout << "nCandidateCut : " << nCandidateCut << endl;
 268     cout << "nCandidateCut2 : " << nCandidateCut2 << endl;
 269     }
 270
 271   // local arrays for input No Cuts
 272   // Both pt < ptMin and pt > ptMin
 273     Float_t* enT       = new Float_t[nCandidate];
 274     Float_t* ptT       = new Float_t[nCandidate];
 275     Float_t* etaT      = new Float_t[nCandidate];
 276     Float_t* phiT      = new Float_t[nCandidate];
 277     Float_t* detaT     = new Float_t[nCandidate];
 278     Float_t* dphiT     = new Float_t[nCandidate];
 279     Float_t* cFlagT    = new Float_t[nCandidate];
 280     Float_t* cClusterT = new Float_t[nCandidate];
 281     Float_t* idT       = new Float_t[nCandidate];
 282     Int_t loop1 = 0;
 283
 284     fJets->SetNinput(nCandidate);
 285
 286     if(fDebug>3){
 287         cout << "nCandidate : " << nCandidate << endl;
 288         //    cout << "nMultCandidate : " << nMultCandidate << endl;
 289     }
 290
 291     Float_t *etCell = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
 292     Float_t *etaCell = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
 293     Float_t *phiCell = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
 294
 295     // Information extracted from fUnitArray
 296     for(Int_t i=0; i<nIn; i++)
 297     {
 298         if(fUnit[i].GetUnitCutFlag()==1){
 299             etCell[i] = fUnit[i].GetUnitEnergy();
 300             if (etCell[i] > 0.0) etCell[i] -= fHeader->GetMinCellEt();
 301             if (etCell[i] < 0.0) etCell[i] = 0.;
 302             etaCell[i] = fUnit[i].GetUnitEta();
 303             phiCell[i] = fUnit[i].GetUnitPhi();
 304         }
 305         else {
 306             etCell[i]  = 0.;
 307             etaCell[i] = fUnit[i].GetUnitEta();
 308             phiCell[i] = fUnit[i].GetUnitPhi();
 309         }
 310
 311         if(fUnit[i].GetUnitEnergy()>0.){
 312             ptT[loop1]   = fUnit[i].GetUnitEnergy();
 313             enT[loop1]   = fUnit[i].GetUnitEnergy();
 314             etaT[loop1]  = fUnit[i].GetUnitEta();
 315             phiT[loop1]  = fUnit[i].GetUnitPhi();
 316             detaT[loop1] = fUnit[i].GetUnitDeta();
 317             dphiT[loop1] = fUnit[i].GetUnitDphi();
 318             cFlagT[loop1]= fUnit[i].GetUnitCutFlag();
 319             idT[loop1]   = fUnit[i].GetUnitID();
 320             loop1++;
 321         }
 322     }
 323
 324
 325     if(fDebug > 40) // For comparison
 326     {
 327         for(Int_t j=0; j<nIn; j++) {
 328             if(etCell[j]>0){
 329                 cout << "etCell[" << j << "] : " << etCell[j] << endl;
 330                 cout << "etaCell[" << j << "] : " << etaCell[j] << endl;
 331                 cout << "phiCell[" << j << "] : " << phiCell[j] << endl;
 332             }
 333         }
 334     }
 335
 336     // Run the algo. Parameters from header
 337     //  Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
 338     Int_t nTot      = nIn;
 339     Float_t minmove = fHeader->GetMinMove();
 340     Float_t maxmove = fHeader->GetMaxMove();
 341     Int_t mode      = fHeader->GetMode();
 342     Float_t precbg  = fHeader->GetPrecBg();
 343     Int_t ierror;
 344
 345     ua1_jet_finder(nIn, nTot, etCell, etaCell, phiCell,
 346                    minmove, maxmove, mode, precbg, ierror);
 347
 348     // sort jets
 349     Int_t * idx  = new Int_t[UA1JETS.njet];
 350     TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
 351
 352     if(fDebug > 20)
 353     {
 354         for(Int_t i = 0; i < UA1JETS.njet; i++)
 355         {
 356             cout << "Number of jets found, UA1JETS.njet : " << UA1JETS.njet << endl;
 357             cout << "UA1JETS.etj : " << UA1JETS.etj << endl;
 358             cout << "idx[" << i << "] : " << idx[i] << endl;
 359             cout << "UA1JETS.etaj[1][" << idx[i] << "] : " << UA1JETS.etaj[1][idx[i]] << endl;
 360             cout << "UA1JETS.phij[1][" << idx[i] << "] : " << UA1JETS.phij[1][idx[i]] << endl;
 361             cout << "UA1JETS.etj[" << idx[i] << "] : " << UA1JETS.etj[idx[i]] << endl;
 362         }
 363     }
 364
 365     // download jet info.
 366     for(Int_t i = 0; i < UA1JETS.njet; i++) {
 367         // reject events outside acceptable eta range
 368         if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
 369             || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
 370         {
 371             continue;
 372         }
 373
 374         Float_t px, py,pz,en,pT; // convert to 4-vector
 375         px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
 376         py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
 377         pz = UA1JETS.etj[idx[i]] /
 378             TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
 379         en = TMath::Sqrt(px * px + py * py + pz * pz);
 380         pT = TMath::Sqrt(px * px + py * py);
 381
 382         fJets->AddJet(px, py, pz, en);
 383
 384     }
 385
 386     // find multiplicities and relationship jet-particle
 387     // find also percentage of pt from pythia
 388
 389     Int_t* injet = new Int_t[nCandidate];
 390     Int_t* sflag = new Int_t[nCandidate];
 391     for (Int_t i = 0; i < nCandidate; i++) {injet[i]= 0;sflag[i]=0;}
 392     Int_t* mult  = new Int_t[fJets->GetNJets()];
 393     Int_t* ncell  = new Int_t[fJets->GetNJets()];
 394     Float_t* percentage  = new Float_t[fJets->GetNJets()];
 395
 396     // Instead of using etaT below, it would be interesting to use the previous fUnitArray object
 397     // With the particle ID, it is possible to to have access to its physical properties and one can,
 398     // for example, set if the corresponding particle is inside or outside the jet with the flag
 399     // kOutJet/kInJet, other possibilities...
 400
 401     for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
 402         Float_t pt_sig = 0.0;
 403         mult[i] = 0;
 404         ncell[i] = UA1JETS.ncellj[i];
 405         for (Int_t j = 0; j < nCandidate; j++) {
 406             Float_t deta = etaT[j] - fJets->GetEta(i);
 407             Float_t dphi = phiT[j] - fJets->GetPhi(i);
 408             if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
 409             if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
 410             Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
 411             if (dr < fHeader->GetRadius() && injet[j] == 0) {
 412                 injet[j] = -(i+1);
 413                 if(cFlagT[j] == 1 &&
 414                    (etaT[j] < fHeader->GetLegoEtaMax()) &&
 415                    (etaT[j] > fHeader->GetLegoEtaMin())) {
 416                     injet[j] = i+1;
 417                     mult[i]++;
 418                     pt_sig+=enT[j];
 419                     sflag[j]=1;
 420                 }
 421             }
 422             if(fDebug>10){
 423                 cout << "mult[" << i << "] : " << mult[i] << endl;
 424                 cout << "ncell[" << i << "] : " << ncell[i] << endl;
 425             }
 426         }
 427         percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
 428             ((Double_t) fJets->GetPt(i));
 429     }
 430
 431     fJets->SetNCells(ncell);
 432     fJets->SetPtFromSignal(percentage);
 433     fJets->SetMultiplicities(mult);
 434     fJets->SetInJet(injet);
 435     fJets->SetEtaIn(etaT);
 436     if(fDebug>10){
 437         for(Int_t i=0; i<nCandidate ; i++){
 438             cout << "phiT[" << i << "] : " << phiT[i] << endl;
 439             cout << "etaT[" << i << "] : " << etaT[i] << endl;
 440         }
 441     }
 442     fJets->SetPhiIn(phiT);
 443     fJets->SetPtIn(enT);
 444     fJets->SetEtAvg(UA1JETS.etavg);
 445     delete etCell;
 446     delete etaCell;
 447     delete phiCell;
 448     delete ncell;
 449     delete cFlagT;
 450     delete cClusterT;
 451     delete enT;
 452     delete ptT;
 453     delete etaT;
 454     delete phiT;
 455     delete detaT;
 456     delete dphiT;
 457     delete injet;
 458     delete idx;
 459     delete mult;
 460     delete percentage;
 461
 462 }
 463
 464
 465 ////////////////////////////////////////////////////////////////////////
 466
 467 void AliUA1JetFinder::Reset()
 468 {
 469   fLego->Reset();
 470   fJets->ClearJets();
 471 }
 472
 473 ////////////////////////////////////////////////////////////////////////
 474
 475 void AliUA1JetFinder::WriteJHeaderToFile()
 476 {
 477   fOut->cd();
 478   fHeader->Write();
 479 }
 480
 481 ////////////////////////////////////////////////////////////////////////
 482
 483 void AliUA1JetFinder::Init()
 484 {
 485   // initializes some variables
 486   Float_t dEta, dPhi;
 487   dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
 488     /((Float_t) fHeader->GetLegoNbinEta());
 489   dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
 490     /((Float_t) fHeader->GetLegoNbinPhi());
 491
 492   UA1CELL.etaCellSize = dEta;
 493   UA1CELL.phiCellSize = dPhi;
 494
 495   // jet parameters
 496   UA1PARA.coneRad = fHeader->GetRadius();
 497   UA1PARA.etSeed  = fHeader->GetEtSeed();
 498   UA1PARA.ejMin   = fHeader->GetMinJetEt();
 499   UA1PARA.etMin   = fHeader->GetMinCellEt();
 500
 501   // book lego
 502   fLego = new
 503     TH2F("legoH","eta-phi",
 504          fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(),
 505          fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(),
 506          fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
 507 }