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 #include "AliAODTrack.h"
  44
  45 #include "fastjet/PseudoJet.hh"
  46 #include "fastjet/ClusterSequenceArea.hh"
  47 #include "fastjet/AreaDefinition.hh"
  48 #include "fastjet/JetDefinition.hh"
  49 // get info on how fastjet was configured
  50 #include "fastjet/config.h"
  51
  52 #ifdef ENABLE_PLUGIN_SISCONE
  53 #include "fastjet/SISConePlugin.hh"
  54 #endif
  55
  56 #include<sstream>  // needed for internal io
  57 #include<vector>
  58 #include <cmath>
  59
  60 using namespace std;
  61
  62
  63 ClassImp(AliFastJetFinder)
  64
  65
  66 //____________________________________________________________________________
  67
  68 AliFastJetFinder::AliFastJetFinder():
  69   AliJetFinder(),
  70   fInputFJ(new AliFastJetInput()),
  71   fJetBkg(new  AliJetBkg())
  72 {
  73   // Constructor
  74 }
  75
  76 //____________________________________________________________________________
  77
  78 AliFastJetFinder::~AliFastJetFinder()
  79 {
  80   // destructor
  81   delete  fInputFJ; fInputFJ = 0;
  82   delete  fJetBkg; fJetBkg = 0;
  83 }
  84
  85 //______________________________________________________________________________
  86 void AliFastJetFinder::FindJets()
  87 {
  88   cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl;
  89   //pick up fastjet header
  90   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  91   Bool_t debug  = header->GetDebug();     // debug option
  92   Bool_t bgMode = header->GetBGMode();    // choose to subtract BG or not
  93
  94   // check if we are reading AOD jets
  95   TRefArray *refs = 0;
  96   Bool_t fromAod = !strcmp(fReader->ClassName(),"AliJetAODReader");
  97   if (fromAod) { refs = fReader->GetReferences(); }
  98
  99   // RUN ALGORITHM
 100   // read input particles -----------------------------
 101
 102   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
 103   if(inputParticles.size()==0){
 104     Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
 105     return;
 106   }
 107
 108   // create an object that represents your choice of jet algorithm, and
 109   // the associated parameters
 110   double rParam = header->GetRparam();
 111   fastjet::Strategy strategy = header->GetStrategy();
 112   fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
 113   fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
 114   fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy);
 115
 116   // create an object that specifies how we to define the area
 117   fastjet::AreaDefinition areaDef;
 118   double ghostEtamax = header->GetGhostEtaMax();
 119   double ghostArea   = header->GetGhostArea();
 120   int    activeAreaRepeats = header->GetActiveAreaRepeats();
 121
 122   // now create the object that holds info about ghosts
 123   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
 124   // and from that get an area definition
 125   fastjet::AreaType areaType = header->GetAreaType();
 126   areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 127
 128
 129
 130
 131
 132   if(bgMode) // BG subtraction!!!!!!!!!! Not entering here
 133     {
 134       //***************************** JETS FINDING AND EXTRACTION
 135       // run the jet clustering with the above jet definition
 136       fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 137
 138       // save a comment in the header
 139
 140       TString comment = "Running FastJet algorithm with the following setup. ";
 141       comment+= "Jet definition: ";
 142       comment+= TString(jetDef.description());
 143       comment+= ". Area definition: ";
 144       comment+= TString(areaDef.description());
 145       comment+= ". Strategy adopted by FastJet: ";
 146       comment+= TString(clust_seq.strategy_string());
 147       header->SetComment(comment);
 148       if(debug){
 149         cout << "--------------------------------------------------------" << endl;
 150         cout << comment << endl;
 151         cout << "--------------------------------------------------------" << endl;
 152       }
 153       //header->PrintParameters();
 154
 155
 156       // extract the inclusive jets with pt > ptmin, sorted by pt
 157       double ptmin = header->GetPtMin();
 158       vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 159
 160       //cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
 161
 162
 163       //subtract background // ===========================================
 164       // set the rapididty , phi range within which to study the background
 165       double rapMax = header->GetRapMax();
 166       double rapMin = header->GetRapMin();
 167       double phiMax = header->GetPhiMax();
 168       double phiMin = header->GetPhiMin();
 169       fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
 170
 171       // subtract background
 172       vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);
 173
 174       // print out
 175       //cout << "Printing inclusive sub jets with pt > "<< ptmin<<" GeV\n";
 176       //cout << "---------------------------------------\n";
 177       //cout << endl;
 178       //printf(" ijet   rap      phi        Pt         area  +-   err\n");
 179
 180       // sort jets into increasing pt
 181       vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
 182       for (size_t j = 0; j < jets.size(); j++) { // loop for jets
 183
 184         double area     = clust_seq.area(jets[j]);
 185         double areaError = clust_seq.area_error(jets[j]);
 186
 187         if(debug) 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);
 188
 189         // go to write AOD  info
 190         AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
 191         //cout << "Printing jet " << endl;
 192         if(debug) aodjet.Print("");
 193         //cout << "Adding jet ... " ;
 194         AddJet(aodjet);
 195         //cout << "added \n" << endl;
 196
 197       }
 198     }
 199
 200
 201
 202
 203
 204
 205   else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
 206
 207     TClonesArray* fUnit = fReader->GetUnitArray(); //Big mmentum array
 208     if(fUnit == 0) { cout << "Could not get the momentum array" << endl; return; }
 209     Int_t         nIn = fUnit->GetEntries();
 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       if(debug) 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       // go to write AOD  info
 246       AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
 247       aodjet.SetEffArea(area,0);
 248       //cout << "Printing jet " << endl;
 249       if(debug) aodjet.Print("");
 250       Int_t count=0;
 251       for (int i=0; i < nCon; i++)
 252         {
 253         fastjet::PseudoJet mPart=constituents[i];
 254         ind[i]=mPart.user_index();
 255         //      cout<<i<<"  index="<<ind[i]<<endl;
 256
 257         //internal loop over all the unit cells
 258         Int_t ipart = 0;
 259         //      count=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                 TRefArray* trackArray = (TRefArray*)uArray->GetUnitTrackRef();
 267                 Int_t tracksInCell = trackArray->GetEntries();
 268                 for(int j=0;j<tracksInCell;j++){
 269                   AliAODTrack * track = (AliAODTrack*)trackArray->At(j);
 270                   aodjet.AddTrack(track);
 271                 }
 272
 273                 count++;
 274               }
 275               ipart++;
 276             }
 277           }
 278       }
 279
 280
 281
 282       AddJet(aodjet);
 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 }