JETAN/DEV/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 /* $Id$ */
  17
  18 //---------------------------------------------------------------------
  19 // FastJet v2.3.4 finder algorithm interface
  20 // Last modification: Neutral cell energy included in the jet reconstruction
  21 //
  22 // Authors: Rafael.Diaz.Valdes@cern.ch
  23 //          Magali.estienne@subatech.in2p3.fr (neutral part + bg subtraction option)
  24 //
  25 // ** 2011 magali.estienne@subatech.in2p3.fr &  alexandre.shabetai@cern.ch
  26 // new implementation of background subtraction
  27 // allowing to subtract bkg using a different algo than the one used for signal jets
  28 //---------------------------------------------------------------------
  29
  30 #include <Riostream.h>
  31
  32 #include "AliFastJetFinder.h"
  33 #include "AliFastJetHeaderV1.h"
  34 #include "AliFastJetInput.h"
  35 #include "AliFastJetBkg.h"
  36 #include "AliAODJetEventBackground.h"
  37 #include "AliAODJet.h"
  38
  39 #include "fastjet/PseudoJet.hh"
  40 #include "fastjet/ClusterSequenceArea.hh"
  41 #include "fastjet/AreaDefinition.hh"
  42 #include "fastjet/JetDefinition.hh"
  43
  44 #include<vector>
  45
  46 using namespace std;
  47
  48 ClassImp(AliFastJetFinder)
  49
  50 ////////////////////////////////////////////////////////////////////////
  51
  52 AliFastJetFinder::AliFastJetFinder():
  53   AliJetFinder(),
  54   fInputFJ(new AliFastJetInput()),
  55   fJetBkg(new  AliFastJetBkg())
  56 {
  57   // Constructor
  58 }
  59
  60 //____________________________________________________________________________
  61 AliFastJetFinder::~AliFastJetFinder()
  62 {
  63   // destructor
  64   delete  fInputFJ;
  65   delete  fJetBkg;
  66
  67 }
  68
  69 //______________________________________________________________________________
  70 void AliFastJetFinder::FindJets()
  71 {
  72   // runs a FASTJET based algo
  73
  74   //pick up fastjet header
  75   AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
  76   Int_t  debug  = header->GetDebug();     // debug option
  77   Bool_t bgMode = header->GetBGMode();    // choose to subtract BG or not
  78   if(debug>0) cout<<"----------in AliFastJetFinder::FindJets() ------------------"<<endl;
  79
  80   // RUN ALGORITHM
  81   // read input particles -----------------------------
  82
  83   vector<fastjet::PseudoJet> inputParticles=fInputFJ->GetInputParticles();
  84   if(inputParticles.size()==0){
  85     if(debug>0) Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
  86     return;
  87   }
  88
  89   // create an object that represents your choice of jet algorithm, and
  90   // the associated parameters
  91   double rParam = header->GetRparam();
  92   double rBkgParam = header->GetRparamBkg();
  93   fastjet::Strategy strategy = header->GetStrategy();
  94   fastjet::RecombinationScheme recombScheme = header->GetRecombScheme();
  95   fastjet::JetAlgorithm algorithm = header->GetAlgorithm();
  96   fastjet::JetDefinition jetDef(algorithm, rParam, recombScheme, strategy);
  97
  98   // create an object that specifies how we to define the area
  99   fastjet::AreaDefinition areaDef;
 100   double ghostEtamax       = header->GetGhostEtaMax();
 101   double ghostArea         = header->GetGhostArea();
 102   int    activeAreaRepeats = header->GetActiveAreaRepeats();
 103
 104   // now create the object that holds info about ghosts
 105   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
 106   // and from that get an area definition
 107   fastjet::AreaType areaType = header->GetAreaType();
 108   areaDef = fastjet::AreaDefinition(areaType,ghostSpec);
 109
 110   //***************************** JETS FINDING
 111   // run the jet clustering with the above jet definition
 112   fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 113
 114   vector<fastjet::PseudoJet> jets;
 115
 116   if(bgMode) // Do BG subtraction directly with the same algorithm (cambridge or kt) for jet signal and background
 117     {
 118       //***************************** JETS FINDING FOR RHO ESTIMATION
 119       // run the jet clustering with the above jet definition
 120       fastjet::JetAlgorithm algorithmBkg = header->GetBGAlgorithm();
 121       fastjet::JetDefinition jetDefBkg(algorithmBkg, rBkgParam, recombScheme, strategy);
 122       fastjet::ClusterSequenceArea clust_seq_bkg(inputParticles, jetDefBkg, areaDef);
 123
 124       // save a comment in the header
 125       TString comment = "Running FastJet algorithm with the following setup. ";
 126       comment+= "Jet definition: ";
 127       comment+= TString(jetDef.description());
 128       comment+= "Jet bckg definition: ";
 129       comment+= TString(jetDefBkg.description());
 130       comment+= ". Area definition: ";
 131       comment+= TString(areaDef.description());
 132       comment+= ". Strategy adopted by FastJet: ";
 133       comment+= TString(clust_seq.strategy_string());
 134       header->SetComment(comment);
 135       if(debug>0){
 136         cout << "--------------------------------------------------------" << endl;
 137         cout << comment << endl;
 138         cout << "--------------------------------------------------------" << endl;
 139       }
 140
 141       // extract the inclusive jets sorted by pt
 142       double ptmin = header->GetPtMin();
 143       vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets();
 144
 145       //subtract background // ===========================================
 146       // set the rapididty , phi range within which to study the background
 147       double rapMax = header->GetRapMax();
 148       double rapMin = header->GetRapMin();
 149       double phiMax = header->GetPhiMax();
 150       double phiMin = header->GetPhiMin();
 151       fastjet::RangeDefinition range(rapMin, rapMax, phiMin, phiMax);
 152
 153       // Extract rho and sigma
 154       Double_t rho = 0.;
 155       Double_t sigma = 0.;
 156       Double_t meanarea = 0.;
 157       Bool_t kUse4VectorArea = header->Use4VectorArea();
 158       vector<fastjet::PseudoJet> bkgJets = clust_seq_bkg.inclusive_jets();
 159       clust_seq_bkg.get_median_rho_and_sigma(bkgJets,range, kUse4VectorArea, rho, sigma, meanarea, false);
 160
 161       // subtract background and extract jets bkg subtracted
 162       vector<fastjet::PseudoJet> subJets = clust_seq.subtracted_jets(rho,ptmin);
 163
 164       // sort jets into increasing pt
 165       jets = sorted_by_pt(subJets);
 166
 167     }
 168
 169   else { // No BG subtraction!!!!!!!! Default header is bgmode=0.
 170
 171     // save a comment in the header
 172     TString comment = "Running FastJet algorithm with the following setup. ";
 173     comment+= "Jet definition: ";
 174     comment+= TString(jetDef.description());
 175     comment+= ". Strategy adopted by FastJet: ";
 176     comment+= TString(clust_seq.strategy_string());
 177     header->SetComment(comment);
 178     if(debug>0){
 179       cout << "--------------------------------------------------------" << endl;
 180       cout << comment << endl;
 181       cout << "--------------------------------------------------------" << endl;
 182     }
 183
 184     // extract the inclusive jets with pt > ptmin, sorted by pt
 185     double ptmin = header->GetPtMin();
 186     vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 187
 188     jets = sorted_by_pt(inclusiveJets);
 189
 190   }
 191
 192   for (size_t j = 0; j < jets.size(); j++) { // loop for jets
 193
 194     double area      = clust_seq.area(jets[j]);
 195     double areaError = clust_seq.area_error(jets[j]);
 196
 197     if(debug>0) printf("Jet found %5d %9.5f %8.5f %10.3f \n",(Int_t)j,jets[j].rap(),jets[j].phi(),jets[j].perp());
 198
 199     vector<fastjet::PseudoJet> constituents = clust_seq.constituents(jets[j]);
 200     int nCon= constituents.size();
 201     TArrayI ind(nCon);
 202
 203     if ((jets[j].eta() > (header->GetJetEtaMax())) ||
 204         (jets[j].eta() < (header->GetJetEtaMin())) ||
 205         (jets[j].phi() > (header->GetJetPhiMax())) ||
 206         (jets[j].phi() < (header->GetJetPhiMin())) ||
 207         (jets[j].perp() < header->GetPtMin())) continue; // acceptance eta range and etmin
 208
 209     // go to write AOD  info
 210     AliAODJet aodjet (jets[j].px(), jets[j].py(), jets[j].pz(), jets[j].E());
 211     aodjet.SetEffArea(area,areaError);
 212     //cout << "Printing jet " << endl;
 213     if(debug>0) aodjet.Print("");
 214
 215     for (int i=0; i < nCon; i++)
 216       {
 217         fastjet::PseudoJet mPart=constituents[i];
 218         ind[i]=mPart.user_index();
 219
 220         // Jet constituents (charged tracks) added to the AliAODJet
 221         AliJetCalTrkEvent* calEvt  = GetCalTrkEvent();
 222         for(Int_t itrack=0; itrack<calEvt->GetNCalTrkTracks(); itrack++)
 223           {
 224             if(itrack==ind[i])
 225               {
 226                 TObject *track = calEvt->GetCalTrkTrack(itrack)->GetTrackObject();
 227                 aodjet.AddTrack(track);
 228               }
 229           }
 230       } // End loop on Constituents
 231
 232     AddJet(aodjet);
 233
 234   } // end loop for jets
 235
 236 }
 237
 238 //____________________________________________________________________________
 239 void AliFastJetFinder::RunTest(const char* datafile)
 240 {
 241   // This simple test run the kt algorithm for an ascii file testdata.dat
 242
 243   // read input particles -----------------------------
 244   vector<fastjet::PseudoJet> inputParticles;
 245   Float_t px,py,pz,en;
 246   ifstream in;
 247   Int_t nlines = 0;
 248   // we assume a file basic.dat in the current directory
 249   // this file has 3 columns of float data
 250   in.open(datafile);
 251   while (1) {
 252     in >> px >> py >> pz >> en;
 253     if (!in.good()) break;
 254     //printf("px=%8f, py=%8f, pz=%8fn",px,py,pz);
 255     nlines++;
 256     inputParticles.push_back(fastjet::PseudoJet(px,py,pz,en));
 257   }
 258   //printf(" found %d pointsn",nlines);
 259   in.close();
 260   //////////////////////////////////////////////////
 261
 262   // create an object that represents your choice of jet algorithm, and
 263   // the associated parameters
 264   double rParam = 1.0;
 265   fastjet::Strategy strategy = fastjet::Best;
 266   fastjet::RecombinationScheme recombScheme = fastjet::BIpt_scheme;
 267   fastjet::JetDefinition jetDef(fastjet::kt_algorithm, rParam, recombScheme, strategy);
 268
 269   // create an object that specifies how we to define the area
 270   fastjet::AreaDefinition areaDef;
 271   double ghostEtamax = 7.0;
 272   double ghostArea    = 0.05;
 273   int    activeAreaRepeats = 1;
 274
 275   // now create the object that holds info about ghosts
 276   fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea);
 277   // and from that get an area definition
 278   areaDef = fastjet::AreaDefinition(fastjet::active_area,ghostSpec);
 279
 280   // run the jet clustering with the above jet definition
 281   fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef);
 282
 283   // tell the user what was done
 284   cout << "--------------------------------------------------------" << endl;
 285   cout << "Jet definition was: " << jetDef.description() << endl;
 286   cout << "Area definition was: " << areaDef.description() << endl;
 287   cout << "Strategy adopted by FastJet was "<< clust_seq.strategy_string()<<endl<<endl;
 288   cout << "--------------------------------------------------------" << endl;
 289
 290   // extract the inclusive jets with pt > 5 GeV, sorted by pt
 291   double ptmin = 5.0;
 292   vector<fastjet::PseudoJet> inclusiveJets = clust_seq.inclusive_jets(ptmin);
 293
 294   cout << "Number of unclustered particles: " << clust_seq.unclustered_particles().size() << endl;
 295
 296   //subtract background // ===========================================
 297   // set the rapididty range within which to study the background
 298   double rapMax = ghostEtamax - rParam;
 299   fastjet::RangeDefinition range(rapMax);
 300   // subtract background
 301   vector<fastjet::PseudoJet> subJets =  clust_seq.subtracted_jets(range,ptmin);
 302
 303   // print them out //================================================
 304   cout << "Printing inclusive jets  after background subtraction \n";
 305   cout << "------------------------------------------------------\n";
 306   // sort jets into increasing pt
 307   vector<fastjet::PseudoJet> jets = sorted_by_pt(subJets);
 308
 309   printf(" ijet   rap      phi        Pt         area  +-   err\n");
 310   for (size_t j = 0; j < jets.size(); j++) {
 311     double area      = clust_seq.area(jets[j]);
 312     double areaError = clust_seq.area_error(jets[j]);
 313     printf("%5d %9.5f %8.5f %10.3f %8.3f +- %6.3f\n",(Int_t)j,jets[j].rap(),
 314            jets[j].phi(),jets[j].perp(), area, areaError);
 315   }
 316   cout << endl;
 317   // ================================================================
 318
 319 }
 320
 321 //____________________________________________________________________________
 322 void AliFastJetFinder::WriteJHeaderToFile() const
 323 {
 324   // Write Jet Header to file
 325   fHeader->Write();
 326
 327 }
 328
 329 //____________________________________________________________________________
 330 Bool_t AliFastJetFinder::ProcessEvent()
 331 {
 332   // Process one event
 333   // Charged only or charged+neutral jets
 334
 335   fInputFJ->SetHeader(fHeader);
 336   fInputFJ->SetCalTrkEvent(GetCalTrkEvent());
 337   fInputFJ->FillInput();
 338
 339   // Find Jets
 340   FindJets();
 341
 342   // Background
 343   if(fAODEvBkg){
 344     AliFastJetHeaderV1 *header = (AliFastJetHeaderV1*)fHeader;
 345
 346     fJetBkg->SetHeader(fHeader);
 347     fJetBkg->SetFastJetInput(fInputFJ);
 348
 349      Int_t count = 0;
 350      if(header->GetBkgFastJetb()){
 351         Double_t sigma1 = 0 , meanarea1= 0, bkg1 = 0;
 352         fJetBkg->BkgFastJetb(bkg1,sigma1,meanarea1);
 353         fAODEvBkg->SetBackground(count,bkg1,sigma1,meanarea1);
 354         count++;
 355      }
 356
 357      if(header->GetBkgFastJetWoHardest()){
 358         Double_t sigma2 = 0 , meanarea2 = 0, bkg2 = 0;
 359         fJetBkg->BkgFastJetWoHardest(bkg2,sigma2,meanarea2);
 360         fAODEvBkg->SetBackground(count,bkg2,sigma2,meanarea2);
 361         count++;
 362      }
 363
 364   }
 365
 366   Reset();
 367   return kTRUE;
 368
 369 }