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