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