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 //---------------------------------------------------------------------
  18 // FastJet v2.3.4 finder algorithm interface
  19 // Last modification: Neutral cell energy included in the jet reconstruction
  20 //
  21 // Authors: Rafael.Diaz.Valdes@cern.ch
  22 //          Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option)
  23 //
  24 //---------------------------------------------------------------------
  25
  26
  27 #include <Riostream.h>
  28 #include <TLorentzVector.h>
  29 #include <TFile.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TArrayF.h>
  33 #include <TClonesArray.h>
  34
  35 #include "AliFastJetFinder.h"
  36 #include "AliFastJetHeaderV1.h"
  37 #include "AliJetReaderHeader.h"
  38 #include "AliJetReader.h"
  39 #include "AliJetUnitArray.h"
  40 #include "AliFastJetInput.h"
  41 #include "AliJetBkg.h"
  42 #include "AliAODJetEventBackground.h"
  43
  44 #include "fastjet/PseudoJet.hh"
  45 #include "fastjet/ClusterSequenceArea.hh"
  46 #include "fastjet/AreaDefinition.hh"
  47 #include "fastjet/JetDefinition.hh"
  48 // get info on how fastjet was configured
  49 #include "fastjet/config.h"
  50
  51 #ifdef ENABLE_PLUGIN_SISCONE
  52 #include "fastjet/SISConePlugin.hh"
  53 #endif
  54
  55 #include<sstream>  // needed for internal io
  56 #include<vector>
  57 #include <cmath>
  58
  59 using namespace std;
  60
  61
  62 ClassImp(AliFastJetFinder)
  63
  64
  65 //____________________________________________________________________________
  66
  67 AliFastJetFinder::AliFastJetFinder():
  68   AliJetFinder(),
  69   fInputFJ(new AliFastJetInput()),
  70   fJetBkg(new  AliJetBkg())
  71 {
  72   // Constructor
  73 }
  74
  75 //____________________________________________________________________________
  76
  77 AliFastJetFinder::~AliFastJetFinder()
  78 {
  79   // destructor
  80   delete  fInputFJ; fInputFJ = 0;
  81   delete  fJetBkg; fJetBkg = 0;
  82 }
  83
  84 //______________________________________________________________________________
  85 void AliFastJetFinder::FindJets()
  86 {
  87   cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl;
  88   //pick up fastjet header
  89   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  90   Bool_t debug  = header->GetDebug();     // debug option
  91   Bool_t bgMode = header->GetBGMode();    // choose to subtract BG or not
  92
  93   // check if we are reading AOD jets
  94   TRefArray *refs = 0;
  95   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
  96   if (fromAod) { refs = fReader->GetReferences(); }
  97
  98   // RUN ALGORITHM
  99   // read input particles -----------------------------
 100
 101   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
 102   if(inputParticles.size()==0){
 103     Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 104     return;
 105   }
 106
 107   // create an object that represents your choice of jet algorithm, and
 108   // the associated parameters
 109   double rParam = header->GetRparam();
 110   fastjet::Strategy strategy = header->GetStrategy();
 111   fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
 112   fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
 113   fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy);
 114
 115   // create an object that specifies how we to define the area
 116   fastjet::AreaDefinition areaDef;
 117   double ghostEtamax = header->GetGhostEtaMax();
 118   double ghostArea   = header->GetGhostArea();
 119   int    activeAreaRepeats = header->GetActiveAreaRepeats();
 120
 121   // now create the object that holds info about ghosts
 122   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
 123   // and from that get an area definition
 124   fastjet::AreaType areaType = header->GetAreaType();
 125   areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 126
 127   if(bgMode) // BG subtraction
 128     {
 129       //***************************** JETS FINDING AND EXTRACTION
 130       // run the jet clustering with the above jet definition
 131       fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 132
 133       // save a comment in the header
 134
 135       TString comment = "Running FastJet algorithm with the following setup. ";
 136       comment+= "Jet definition: ";
 137       comment+= TString(jetDef.description());
 138       comment+= ". Area definition: ";
 139       comment+= TString(areaDef.description());
 140       comment+= ". Strategy adopted by FastJet: ";
 141       comment+= TString(clust_seq.strategy_string());
 142       header->SetComment(comment);
 143       if(debug){
 144         cout << "--------------------------------------------------------" << endl;
 145         cout << comment << endl;
 146         cout << "--------------------------------------------------------" << endl;
 147       }
 148       //header->PrintParameters();
 149
 150
 151       // extract the inclusive jets with pt > ptmin, sorted by pt
 152       double ptmin = header->GetPtMin();
 153       vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 154
 155       //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
 156
 157
 158       //subtract background // ===========================================
 159       // set the rapididty , phi range within which to study the background
 160       double rapMax = header->GetRapMax();
 161       double rapMin = header->GetRapMin();
 162       double phiMax = header->GetPhiMax();
 163       double phiMin = header->GetPhiMin();
 164       fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
 165
 166       // subtract background
 167       vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);
 168
 169       // print out
 170       //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n";
 171       //cout << "---------------------------------------\n";
 172       //cout << endl;
 173       //printf(" ijet   rap      phi        Pt         area  +-   err\n");
 174
 175       // sort jets into increasing pt
 176       vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
 177       for (size_t j = 0; j < jets.size(); j++) { // loop for jets
 178
 179         double area     = clust_seq.area(jets[j]);
 180         double areaError = clust_seq.area_error(jets[j]);
 181
 182         printf("Jet found %5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n", (Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, areaError);
 183
 184         // go to write AOD  info
 185         AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
 186         //cout << "Printing jet " << endl;
 187         if(debug) aodjet.Print("");
 188         //cout << "Adding jet ... " ;
 189         AddJet(aodjet);
 190         //cout << "added \n" << endl;
 191
 192       }
 193     }
 194   else { // No BG subtraction
 195
 196     TClonesArray* fUnit = fReader->GetUnitArray();
 197     if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
 198     Int_t         nIn = fUnit->GetEntries();
 199     cout<<"===== check Unit Array in AliFastJetFinder ========="<<endl;
 200     Int_t ppp=0;
 201     for(Int_t ii=0; ii<nIn; ii++)
 202       {
 203         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(ii);
 204         if(uArray->GetUnitEnergy()>0.){
 205           Float_t eta   = uArray->GetUnitEta();
 206           Float_t phi   = uArray->GetUnitPhi();
 207           cout<<"ipart = "<<ppp<<" eta="<<eta<<"  phi="<<phi<<endl;
 208           ppp++;
 209         }
 210       }
 211
 212     //fastjet::ClusterSequence clust_seq(inputParticles, jetDef);
 213     fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 214
 215     // save a comment in the header
 216
 217     TString comment = "Running FastJet algorithm with the following setup. ";
 218     comment+= "Jet definition: ";
 219     comment+= TString(jetDef.description());
 220     comment+= ". Strategy adopted by FastJet: ";
 221     comment+= TString(clust_seq.strategy_string());
 222     header->SetComment(comment);
 223     if(debug){
 224       cout << "--------------------------------------------------------" << endl;
 225       cout << comment << endl;
 226       cout << "--------------------------------------------------------" << endl;
 227     }
 228     //header->PrintParameters();
 229
 230       // extract the inclusive jets with pt > ptmin, sorted by pt
 231     double ptmin = header->GetPtMin();
 232     vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 233
 234     //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
 235
 236     vector<fastjet::PseudoJet> jets = sorted_by_pt(inclusiveJets); // Added by me
 237     for (size_t j = 0; j < jets.size(); j++) { // loop for jets     // Added by me
 238
 239       printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp());
 240
 241       vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
 242       int nCon= constituents.size();
 243       TArrayI ind(nCon);
 244       Double_t area=clust_seq.area(jets[j]);
 245       cout<<"area = "<<area<<endl;
 246       // go to write AOD  info
 247       AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
 248       aodjet.SetEffArea(area,0);
 249       //cout << "Printing jet " << endl;
 250       if(debug) aodjet.Print("");
 251       // cout << "Adding jet ... " <<j<<endl;
 252       for (int i=0; i < nCon; i++)
 253         {
 254         fastjet::PseudoJet mPart=constituents[i];
 255         ind[i]=mPart.user_index();
 256         //cout<<i<<"  index="<<ind[i]<<endl;
 257
 258         //internal oop over all the unit cells
 259         Int_t ipart = 0;
 260         for(Int_t ii=0; ii<nIn; ii++)
 261           {
 262             AliJetUnitArray *uArray = (AliJetUnitArray*)fUnit->At(ii);
 263             if(uArray->GetUnitEnergy()>0.){
 264               uArray->SetUnitTrackID(ipart);//used to have the index available in AliJEtBkg
 265               if(ipart==ind[i]){
 266                 aodjet.AddTrack(uArray);
 267               }
 268               ipart++;
 269             }
 270           }
 271       }
 272
 273       AddJet(aodjet);
 274       //cout << "added \n" << endl;
 275
 276
 277       ///----> set in the aod the reference to the unit cells
 278       //  in FastJetFinder: 1) loop over the unit array. 2) select those unit cells belonging to the jet (via user_index). 3) use AliAODJet->AddTrack(unitRef)
 279       //  in AliJetBkg: 1) loop over the unit arrays --> get ind of the unit cell 2) internal loop on jet unit cells; 3) check if i_cell = UID of the trackRefs of the AODJet
 280       // should work hopefully
 281
 282
 283
 284     } // end loop for jets
 285   }
 286
 287 }
 288
 289 //____________________________________________________________________________
 290 void AliFastJetFinder::RunTest(const char* datafile)
 291
 292 {
 293
 294    // This simple test run the kt algorithm for an ascii file testdata.dat
 295    // read input particles -----------------------------
 296   vector<fastjet::PseudoJet> inputParticles;
 297   Float_t px,py,pz,en;
 298   ifstream in;
 299   Int_t nlines = 0;
 300   // we assume a file basic.dat in the current directory
 301   // this file has 3 columns of float data
 302   in.open(datafile);
 303   while (1) {
 304       in >> px >> py >> pz >> en;
 305       if (!in.good()) break;
 306       //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz);
 307       nlines++;
 308       inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en));
 309    }
 310    //printf(" found %d pointsn",nlines);
 311    in.close();
 312    //////////////////////////////////////////////////
 313
 314   // create an object that represents your choice of jet algorithm, and
 315   // the associated parameters
 316   double rParam = 1.0;
 317   fastjet::Strategy strategy = fastjet::Best;
 318   fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme;
 319   fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy);
 320
 321
 322
 323   // create an object that specifies how we to define the area
 324   fastjet::AreaDefinition areaDef;
 325   double ghostEtamax = 7.0;
 326   double ghostArea    = 0.05;
 327   int    activeAreaRepeats = 1;
 328
 329
 330   // now create the object that holds info about ghosts
 331   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
 332   // and from that get an area definition
 333   areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec);
 334
 335
 336   // run the jet clustering with the above jet definition
 337   fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 338
 339
 340   // tell the user what was done
 341   cout << "--------------------------------------------------------" << endl;
 342   cout << "Jet definition was: " << jetDef.description() << endl;
 343   cout << "Area definition was: " << areaDef.description() << endl;
 344   cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl;
 345   cout << "--------------------------------------------------------" << endl;
 346
 347
 348   // extract the inclusive jets with pt > 5 GeV, sorted by pt
 349   double ptmin = 5.0;
 350   vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 351
 352   cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
 353
 354
 355   //subtract background // ===========================================
 356   // set the rapididty range within which to study the background
 357   double rapMax = ghostEtamax - rParam;
 358   fastjet::RangeDefinition range(rapMax);
 359   // subtract background
 360   vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);
 361
 362   // print them out //================================================
 363   cout << "Printing inclusive jets  after background subtraction \n";
 364   cout << "------------------------------------------------------\n";
 365   // sort jets into increasing pt
 366   vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
 367
 368   printf(" ijet   rap      phi        Pt         area  +-   err\n");
 369   for (size_t j = 0; j < jets.size(); j++) {
 370
 371     double area     = clust_seq.area(jets[j]);
 372     double areaError = clust_seq.area_error(jets[j]);
 373
 374     printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(),
 375            jets[j].phi(),jets[j].perp(), area, areaError);
 376   }
 377   cout << endl;
 378   // ================================================================
 379
 380
 381
 382 }
 383
 384 //____________________________________________________________________________
 385
 386 void AliFastJetFinder::WriteJHeaderToFile()
 387 {
 388   fHeader->Write();
 389 }
 390
 391 //____________________________________________________________________________
 392
 393 Float_t  AliFastJetFinder::EtaToTheta(Float_t arg)
 394 {
 395   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
 396   return 2.*atan(exp(-arg));
 397
 398
 399 }
 400
 401 //____________________________________________________________________________
 402
 403 void AliFastJetFinder::InitTask(TChain *tree)
 404 {
 405
 406   printf("Fast jet finder initialization ******************");
 407   fReader->CreateTasks(tree);
 408
 409 }
 410
 411
 412 Bool_t AliFastJetFinder::ProcessEvent()
 413 {
 414   //
 415   // Process one event
 416   // from meomntum array
 417
 418   Bool_t ok = fReader->FillMomentumArray();
 419
 420   if (!ok) return kFALSE;
 421   fInputFJ->SetHeader(fHeader);
 422   fInputFJ->SetReader(fReader);
 423   fInputFJ->FillInput();
 424   // Jets
 425   FindJets();
 426
 427   fJetBkg->SetHeader(fHeader);
 428   fJetBkg->SetReader(fReader);
 429   /*
 430   fJetBkg->SetFastJetInput(fInputFJ);
 431   Double_t bkg1=fJetBkg->BkgFastJet();
 432   Double_t bkg2=fJetBkg->BkgChargedFastJet();
 433   Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets);
 434   Double_t bkg4=fJetBkg->BkgRemoveJetLeading(fAODjets);
 435
 436   fAODEvBkg->SetBackground(0,bkg1);
 437   fAODEvBkg->SetBackground(1,bkg2);
 438   fAODEvBkg->SetBackground(2,bkg3);
 439   fAODEvBkg->SetBackground(3,bkg4);
 440   */
 441   Reset();
 442   return kTRUE;
 443
 444 }
 445
 446 Bool_t AliFastJetFinder::ProcessEvent2()
 447 {
 448   //
 449   // Process one event
 450   // Charged only or charged+neutral jets
 451   //
 452
 453   TRefArray* ref = new TRefArray();
 454   Bool_t procid = kFALSE;
 455   Bool_t ok = fReader->ExecTasks(procid,ref);
 456
 457   // Delete reference pointer
 458   if (!ok) {delete ref; return kFALSE;}
 459
 460   // Leading particles
 461   fInputFJ->SetHeader(fHeader);
 462   fInputFJ->SetReader(fReader);
 463   fInputFJ->FillInput();
 464
 465   // Jets
 466   FindJets();
 467
 468   fJetBkg->SetHeader(fHeader);
 469   fJetBkg->SetReader(fReader);
 470   fJetBkg->SetFastJetInput(fInputFJ);
 471   Double_t bkg1=fJetBkg->BkgFastJet();
 472   Double_t bkg2=fJetBkg->BkgChargedFastJet();
 473   Double_t bkg3=fJetBkg->BkgFastJetCone(fAODjets);
 474   Double_t bkg4=fJetBkg->BkgRemoveJetLeading(fAODjets);
 475
 476   fAODEvBkg->SetBackground(0,bkg1);
 477   fAODEvBkg->SetBackground(1,bkg2);
 478   fAODEvBkg->SetBackground(2,bkg3);
 479   fAODEvBkg->SetBackground(3,bkg4);
 480
 481   Int_t nEntRef    = ref->GetEntries();
 482
 483   for(Int_t i=0; i<nEntRef; i++)
 484     {
 485       // Reset the UnitArray content which were referenced
 486       ((AliJetUnitArray*)ref->At(i))->SetUnitTrackID(0);
 487       ((AliJetUnitArray*)ref->At(i))->SetUnitEnergy(0.);
 488       ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag(kPtSmaller);
 489       ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag2(kPtSmaller);
 490       ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlag(kBad);
 491       ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlagC(kTRUE,kBad);
 492       ((AliJetUnitArray*)ref->At(i))->SetUnitDetectorFlag(kTpc);
 493       ((AliJetUnitArray*)ref->At(i))->SetUnitFlag(kOutJet);
 494       ((AliJetUnitArray*)ref->At(i))->ClearUnitTrackRef();
 495
 496       // Reset process ID
 497       AliJetUnitArray* uA = (AliJetUnitArray*)ref->At(i);
 498       uA->ResetBit(kIsReferenced);
 499       uA->SetUniqueID(0);
 500     }
 501
 502   // Delete the reference pointer
 503   ref->Delete();
 504   delete ref;
 505
 506   Reset();
 507
 508   return kTRUE;
 509 }